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ABSTRACT 

The detection of large quantities of dust in z ~ 6 quasars by infrared and radio surv eys presents puzzles for the 
formation and evolution of dust in these early systems. Previously dLi et a l. 2007), we showed that luminous 
quasars at z > 6 can form through hierarchical mergers of gas-rich galaxies, and that these systems are expected 
to evolve from starburst through quasar phases. Here, we calculate the dust properties of simulated quasars and 
their progenitors using a three-dimensional Monte Carlo radiative transfer code, ART 2 - All-wavelength Radia- 
tiv e Transfer with Adaptive R efinement Tree. ART 2 incorporates the radiative equilibrium algorithm developed 
bv lBiorkman & Wood! d2001l) which treats dust emission self-consistently, an adaptive grid method which can 
efficiently cover a large dynamic range in spatial scales and can capture inhomogeneous density distributions, 
a multiphase model of the interstellar medium which accounts for the observed scaling relations of molecular 
clouds, and a supernova-origin model for dust which can explain the existenc e of dust in cosmologically young, 
starbursting quasars. By applying ART 2 to the hydrodynamic simulations of iLi et al] d2007l) . we reproduce the 
observed spectral energy distribution (SED) and inferred dust properties of SDSS J 1 148+5251, the most distant 
quasar detected in the Sloan survey. We find that the dust and infrared emission are closely associated with 
the formation and evolution of the quasar host. As the system evolves from a starburst to a quasar, the SED 
changes from being dominated by a cold dust bump (peaking at ~ 50 /im) to one that includes a prominent 
hot dust component (peaking at ~ 3 /im), and the galaxy evolves from a cold to a warm ultraluminous infrared 
galaxy (ULIRG) owing to heating and feedback from star formation and the active galactic nucleus (AGN). 
Furthermore, the AGN activity has significant implications for the interpretation of observable aspects of the 
hosts. The hottest dust (T > 10 3 K) is most noticeable only during the peak quasar activity, and correlates with 
the near-IR flux. However, we find no correlation between the star formation rate and far-IR luminosity during 
this phase owing to strong AGN contamination. Our results suggest that vigorous star formation in merging 
progenitors is necessary to reproduce the observed dust properties of z ~ 6 quasars, supporting a merger-driven 
origin for luminous quasars at high redshifts and the starburst-to-quasar evolutionary hypothesis. 
Subject headings: quasar: formation — quasar: evolution — quasar: high redshift — galaxies: starbursts 

— infrared: galaxies — radiative transfer — interstellar medium — dust, extinction — 

individual: SDSS Jl 148+5251 



1. INTRODUCTION 

High-redshift quasars are important for understanding the 
formation and evolution of galaxies and supermassive black 
holes (SMBHs) in the early Universe. In the past few years, 
nearly two dozen luminous quasars have been discovered 
by the Sloan Digital Sky Survey (SDSS; lYorket al.ll2000h 
and the Canada-F rance High-z Quasar Survey (CFHQS; 
IWillott etal] 12007) at z ~ 6, corresponding to a ti me when 
the Universe was le ss than a billion years o ld dFan et al] 
l200l 120041 l20063hh . As summarized by iFanl (l2~006h . 
these quasars are rare (~ 10~ 9 Mpc~ 3 comoving); believed 
to be powered by SMBHs with masses ~ 10 9 M Q (e.g., 
IWillott et al] 120031: iBarth et all 120031) : reside near the end 
of the epoch of reio nization, as indicated b y Gunn-Peterson 
absorption troughs dGunn & Petersonl fl965) in their spectra; 
and have similar spectral energy distributions (SEDs) and 
comparable metallicity to lower-redshift counterparts (e.g. , 
Elvis et all 1 19941 iGlikman et all 120061: iRichards et al] 120061; 



HopkinsetafJl2007d), implying the early presence of "ma- 
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ture" quasars and the formation of their hosts in rapid star- 
bursts at even higher redshifts (z > 10). 

Follow-up, multi-wavelength observations ha ve been car- 
ried out for these z ~ 6 qua s ars, from X-ray (e.gjBrandt et all 
20021 ; IStrateva et alj 120051; IVignali et al] 120051 ISteffen et al] 
20061 [ Shemmer et al]|2005l 120061). to optical /infrared (e.g. , 
Barthet al. 2003: Pe ntericci et al. 2003; Freudling et alJ2003t 
White et al.ll2 005l: IWillott et al.ll200H iHines et alJl2006l). and 
radio wavelengths (e.g.. [Carilli et all 1200 lj; iBertoldi et al 



IWalter et all 120031: ICarilli'etall 2004; Wang et al. 
A noteworthy result of these studies is their detec- 
tion of dust in these high-redshift obj ects. In particular , 
deep infrared and radio surveys (e.g., iRobson et al.l 120041; 
IBertoldi et alJl2003at ICarilli et ail 120041: ICharmandaris et all 
12004 iBeelen et al.l l2006h IHines et alj 120061) have revealed a 
large amount of cold dust in SDSS Jl 148+5251 (hereafter 
Jl 148+5251), the most distan t Sloan quasar detected at red- 
shift z = 6.42 (iFan et alJ |2003). The dust mass is estimated to 
be ~ 1 -7 x 10 8 M Q . The detection of dust is also reported 
in the first four CFHQS quasars at z > 6, includ ing the new 
recor d holder CFHQS J2329 -0301 at z = 6.43 (IWillott et alJ 
l2007h . iMaiolino et af] (2004) argue, from the observed dust 
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extinction curve of SDSS J1048+4637 at z = 6.2, that the dust 
in these high-z syste ms may be produced by supernovae. 

I Jiang et ail J2006) have performed a comprehensive study 
of thirteen z ~ 6 quasars by combining new Spitzer observa- 
tions with existing multi-band data. It appears that nearly all 
13 of these quasars exhibit prominent infrared bumps around 
both A ~ 3fim and A ~ 50^m. In a dusty galaxy, most of 
the radiation from newly formed stars or an AGN is ab- 
sorbed by dust and re-emitted at infrared wavelengths. The 
SEDs of star-forming galaxies typically peak at 50 - 80 fim 
(rest frame) and can be approximated by a modified black- 
body spectrum with dust temperatures below 100 K (e.g., 
Dunne & Eales 2001), while in quasars, some dust can be 
heate d directly by the AG N to temperatures up to 1200 K 
(e.g., [Glikman et al. 2006) and dominate t he near- to mid-IR 
emission. Therefore, the observations by Jiang et alj (12006b 
indicate the presence of large amounts of both hot and cold 
dust in these high-z systems. 

However, the nature of the dust, and its formation and evo- 
lution in the context of the quasar host, are unclear. For ex- 
ample, the dust distribution is unknown. It has been suggested 
that the hot dust lies in the central re gions and is heated by an 
AGN to produce near-IR emission (|Rieke & Lebofskvlll981b 
iPollettaet all 120001: lHaas et al.1 120031) . while the warm and 
cold dust can extend to a few kpc and dominate at mid-IR and 
far-IR wavelengths (jPolletta et al.ll2000tlNenkova et al.ll2002t 
Siebenmorgen et al. 2005). However, the relative importance 
of heating by stars and AGN activity is uncertain. This is es- 
sential to an accurate determination of the star formation rate 
(SFR), but cannot be inferred simply from observed SEDs. 
Currently, the most common method used to derive a star 
formation rate is to assume that most of the FIR luminosity 
comes from young stars. However, if the FIR luminosity is 
mainly contributed by an AGN, then the SFR would be sub- 
stantially reduced. Finally, it is not known how the SED and 
dust content of a quasar and its host evolve with time. 

In order to address these questions, we must combine a 
model for the formation of a quasar with radiative transfer cal- 
culations that treat dust emission self-consistently, to follow 
the physical properties, environment, and evolution of quasar 
hosts and their dust conte nt. 

Earlier dLi et al.ll2007l) . we developed a quasar formation 
model which self-consistently accounts for black hole growth, 
star formation, quasar activity, and host spheroid formation in 
the context of hierarchical structure formation. We employed 
a set of multi-scale simulations that included large-scale cos- 
mological simulations and galaxy mergers on galactic scales, 
together with a self-regulated model for black hole growth, 
to produce a luminous quasar at z ~ 6.5, which has a black 
hole mass (~ 2 x 10 9 M^) and a number of properties similar 
to Jl 148+5251 (|FmieLaljl2QQl. In our scenario, luminous, 
high-redshift quasars form in massive halos that originate 
from rare density peaks in the standard ACDM cosmology, 
and they grow through hierarchical mergers of gas-rich galax- 
ies. Gravitational torques excited in these mergers trigger 
large inflows of g as and produce strong shocks that result in 
intense starb ursts (Hernquist 1989; Barnes & Hernauis ti 19911 
19921 1 19961 |Hernc}uist & Mihos 1995; Mihos &Hemc|uisd 



1996) and fuel rapid black hole accretion (Di Matte o et al.l 
20051 ISpringel et all l2005al) . Moreover, feedback from 



the accreting black hole disperses the obsc uring material, 
briefly yielding an optically visible quasar (Hopki ns et al.1 
l2005dfblla l 120061) . and regulating the M bh~c correlation be- 
tween the SMBHs and host galaxies (Di Matteo et af] 120051; 



iRobertson et al.l 120061) . In this picture, quasars are, there- 
fore, descendents of starburst galaxies (ISanders et al.l IT988b 
iNorman & Scovilldl 1 988b IScovilldl2003h . 

To make detailed contact with observations of z ~ 6 
quasars, it is necessary to theoretically predict the SEDs of 
these systems and how they evolve with time. Solving ra- 
diation transport in dusty, starbursting quasars is, however, 
difficult owing to the non-locality of the sources, the opacity 
and multiphase character of the interstellar medium, and the 
complex morphology of these systems. Over the last several 
decades, a variety of numerical techniques have been devel- 
oped to solve the radiative transfer problem with different lev- 
els of approximation and in multi-dimensions. Based on the 
specific algorithms employed, the approaches can be classi- 
fied into two general categories: fin ite-differe nce and Monte 
Carlo methods (e.g., iJonssonl 12006b see also iPascucci et all 
|2004| who describe the codes as "grid-based" or "particle- 
based" by analogy to hydrodynamic solvers). 

Finite-difference codes solve the equations of radiative 
transfer (RT) iteratively using finite convergence criteria. 
They either sol ve the moment equat i ons fo r RT, originally 
formulated by Hummer & Rybick^ d!971[) for spherical 
geometry wit h a ce ntral po int so urce (e.g.. IScoville & Kwanl 
1976b iLeund [19761 lYorSfl980l IWolfire & Cassinellilll986b 
Men'shchikov & Henningl 1 1997b iDullemond & Turollal 
2000), or employ ray-tracing methods on a discrete 
spatial grid for comple x density configurations (e.g., 
Rowan-Robinsonl 1 1980b [E fstathiou & Rowan-Robinson 
19901 ll991b~Tsteinacker"etal.l 12003b iFolini et al.l L2003T 



Steinacker et al.1 120061) . This technique provides full error 



and 



control but can be time-consuming 

Monte Carlo method s sample a i 
proba bil istically (e.g., IWittl 1 1977b 
1983 1 IWhitnev & Hartmannl 11992b 
Code &Whitnevl 1 1995b lLopezetaF 
Wolf etal.1 11999b iBianchi et aT~ 



propagate p hotons 
iLefevre et all 1 19821 



IWitt et all 
11995b iLucvl 
2000b iHarriesI 



1992; 
1999] 
200C 



Biorkman & Woodl l200Tb IWhitnev et all l2003bb IJonssonl 
2006: lPinte et al.ll2006l) . The Monte Carlo technique is more 
flexible than the finite-difference one because it tracks the 
scattering, absorption and re-emission of photons (or photon 
packets) in detail, and can handle arbitrary geometries, but at 
the cost of computational expense to reduce Poisson noise. 
However, advances in computing technology and algorithms 
have made high-accuracy Monte Carlo RT calculations 

feasible and po pular. 

In particular, Biorkman & Wood (200H) have developed a 
Monte Carlo code to handle radiative equilibrium and tem- 
perature corrections, which calculates dust emission self- 
consistently. It conserves the total photon energy, corrects 
the frequency distribution of re-emitted photons, and requires 
no iteration as the dust opacity is assumed to be indepen- 
dent of temperature. This code has been used in a num- 
ber of applications, including protostars that have a disk 
and an e nvelope with a single heating source in the cen- 
ter (e.g. , IWhitnev & Hartmannl 1 19921 1 1 1993b IWhitnev et aL 



2003a b, 2004); circumstellar envelopes (e.g., Wood et al 
1 996alfbL 119981): pro toplanetary systems (e.g., IWood et al. 
2002b iRice et ali 120031) : and galaxies that include a bulge 
and a disk with multiple heating sources from these tw o 
populations (e.g., IWood & Jonell 1997b IWood & Loebll2000l) . 
This code has also been able to generate optical-far-IR 
SEDs that re produce those of a samp le of 21 X-ray se- 
lected AGN s (Kuraszkie wicz et ai1l2003l) . and the version of 
IWhitnev et al.l (l2003bl) has been applied to simulations of 
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galaxy m ergers with black hole fe edback to study the local 
UL IRGs jC hakrabarti et al. 2007a) and submillimeter galax- 
ies dChakrabarti et al.ll2007bl) ~ 

In order to produce accurate S EDs of quasars and their hosts 
formed by galaxy mergers, as in lLi et afl (120071) . the RT code 
must satisfy the following requirements: (1) to be able to han- 
dle arbitrary geometries and distributed heating sources; (2) 
resolve the large dynamic ranges in spatial scales and densi- 
ties in the merger simulations; (3) incorporate a multiphase 
description of the interstellar medium (ISM); and (4) employ 
a dust model appropriate for a young starburst system. 

In many earlier applications, radial logarithmic or uni- 
formly spaced meshes were used to generate density grids. 
However, in merging galaxies, the ISM is clumpy and irreg- 
ular owing to shocks and tidal features produced during the 
interactions, and has multiple density centers, making loga- 
rithmic algorithms inefficient. In such a situation, an adap- 
tive grid approach appears ideal for resolving localized high- 
density regions, while still covering the large volumes of 
merging syste ms (Wolf et al. 1999: Kurosawa & Hillierll200 lh 
lHarries et aT]|2004t lJonssonll2006h 

In addition, dust models commonly used in previous 
works are based on observed extinction curves for the Milky 
Way (e.g., IWeingartner & Drainel 120011; ICalzetti et alj|1994t 
iKim et al J [1994 iMathis et all 1 1 9771: ICalzetti et all 120001) . in 
which the dust is assumed to be prod uced mainly by old, 
low-mass stars with ages > 1 Gyr (iMafhisI fl9 90: Whittet 
120031: iDwekl 120051) . However, a large amount of dust 
(~ 1 - 7 x 10 8 M(7)) is de tected in the z ~ 6.42 quasar 
host dBertoldi et alJ l2003al) at a time when the Universe 
was only ~ 850 Myr old. It has been suggested by ob- 
serva t ions (e.g.. |Maiolino et alJ |2004 120061: iMoselev et alJ 
1989tlDunne et al.ll2003HMorgan et atll2003USugerman et all 
20061) and the oretical studies ("e.g.. | Todini & Ferraral 120011: 
NozawaetaT] 120031: fSchneider et all 12004 iNozawaetafl 
20071: iBianchi & Schneider! 120071) that supernovae can pro- 
vide fast and efficient dust formation in the early Universe, 
motivating other choices for the dust model in the RT calcu- 
lations. 

Finally, in a multip hase description of the ISM 
dMcKee & Ostrikerl 1 19771), as ad o pted in our simula - 
tions dSpringel & Hernquistl l2003allbl: iHopkins etail 120061) . 
the "hot-phase" (diffuse) and "cold-phase" (dense molecular 
and HI core) components co-exist under pressure equilibrium 
but have different mass fractions and volume filling factors 
(i.e., the hot-phase gas is < 10% in mass but > 99% in 
volume), so both phases contribute to the dust extinction and 
should therefore be included in the RT calculations. 

We have refined th e Mont e Car lo RT code d e veloped 
by iBjorkman & Wood! 42001) and IWhitney et afl (l2003bh 
by imp lementing : an adaptive grid algorithm similar to 
that of iJonssonl d2006l) for the density field and the ar- 
bitrarily distributed sources ; a multiphase ISM model 
( Springel & Hernquist 2003a) which accounts for observed 
scaling relations of molecular clouds for the dust distribution; 
and a supernova-origin d ust model for the opacit y using the 
grain size distribution of iTodini & Ferraral (1200 ll) . We refer 
to our new code as ART 2 (All-wavelength Radiative Transfer 
with Adaptive Refinement Tree). ART 2 is capable of produc- 
ing SEDs and images in a wide range of wavelengths from 
X-ray to millimeter. In the present paper we focus on the 
dust properties from optical to submillimeter bands. As we 
show in what follows, ART 2 reproduces the spectrum of a 
single galaxy with bulge and disk calculated with the orig- 



inal code. More important, it captures the inhomogeneous 
density distribution in galaxy mergers and reproduces the ob- 
served SED a nd dust properti es of Jl 148+5251 based on the 
simulations of iLi et all d2007l) . Therefore, ART 2 can be used 
to predict multi-wavelength properties of quasar systems and 
their galaxy progenitors. 

This paper is organized as follows. In § 2, we describe our 
computational methods a nd mo dels, including the multi-scale 
simulations of ILi et a l. (2007) for z ~ 6 quasar formation, 
and our implementat ion of ART 2 which incorp orates radia- 
tive equilibrium as in Bio rkman & Wood! d2001l) . an adaptive 
grid, a multiphase ISM, and a supernova-origin dust model. 
In § 3, we present the multi-wavelength SED from optical to 
submillimeter of a simulated quasar at z ~ 6.5. The dust dis- 
tribution and properties of the quasar are discussed in § 4 and 
we describe their evolution in § 5. We consider the robustness 
of our results in § 6 and summarize in § 7. 

2. METHODOLOGY 

In order to capture the physical processes underlying the 
formation and evolution of quasars in the early Universe, and 
to compare their multi-wavelength properties to observations, 
we combine quasar formation and radiative transfer calcu- 
lations. We first perform a set of novel mul ti-scale simula - 
tions that yield a luminous quasar at z ~ 6.5 dLi et al.l l2007). 
We then apply the 3-D, Monte Carlo radiative transfer code, 
ART 2 , to the outputs of the hydrodynamic galaxy mergers 
simulations to calculate the SED of the system. The models 
and simulation s of quasar formation are described in detail in 
ILi et alJ d2007l) . Here, we briefly summarize the modeling of 
quasar formation, and the specifications of ART 2 . 

2.1. Formation Model of z ~ 6 Quasars 
2.1.1. Multi-scale Simulations 

The z ~ 6 quasars are rare (space density ~ 1 Gpc~ 3 comov- 
ing), and appear to be powered by supermassive black holes 
of mass ~ 10 9 M Q . Therefore, simulations of high-redshift 
quasar formation must consider a large cosmological volume 
to accommodate the low abundance of this population, have a 
large dynamic range to follow the hierarchical build-up of the 
quasar hosts, and include realistic prescriptions for star for- 
mation, black hole growth, and associa ted feedback mecha- 
nisms. The multi-scale simulations in lLi et ail d2007l) include 
both N-body cosmological calculations in a volume of 3 Gpc 3 
to account for the low number density of quasars at z ~ 6, and 
hydrodynamical simulations of individual galaxy mergers on 
galactic scales to resolve gas-dynamics, star formation, and 
black hole growth. 

First, we perform a coarse dark matter-only simulation in 
a volume of 3 Gpc 3 . The largest halo at z = 0, within which 
early, l uminous quasars are thought to reside ( Sprin gel et alJ 
2005b), is then selected for resimulation at higher resolu- 
tion. The evolution of this halo and its environ ment is re- 
simulated using a multi-grid zoom-in technique dGao et alJ 
120051) that provides much higher mass and spatial resolu- 
tion for the halo of interest, while following the surround- 
ing large-scale structure at lower resolution. The merging 
history of the largest halo at z ~ 6, which has then reached 
a mass of ~ 7.7 x 10 12 M Q through seven major (mass ra- 
tio < 5:1) mergers between redshifts 14.4 and 6.5, is ex- 
tracted. These major mergers are again re-simula ted hydro- 
dynam ically using galaxy models which include a Hernquist 
( 1990) halo and an exponential disk scaled appropriately for 
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redshift dRobertson et al.ll2.Q06b . and adjusted to account for 
mass accretion through minor mergers. Each of these eight 
galaxy progenitors has a black hole seed assumed to orig- 
inate from the remnant s of the first stars (|Abel et al.l 120021; 
Bromm & Larsonl l2004t fern & McKeel 12004 lYoshida et al. 
2001 l2006t iGao et al.ll2007h . 

The simulations were performed using the parallel, N- 
body /Smo othed P a rticle Hydrodynamics (SPH) code GAD- 
GET2 (Springel 2005), which conserves energy and 
entrop y using the variational pr inciple formulation of 
SPH (Springel & Hernquist 2002), and which incorpo- 
rates a sub-resolution model of a multiphase interstellar 
medium to describe star form ation and supernova feedback 
(Springel & Hernquist 2003a). Star formation is modeled fol- 
lowing the Schmidt-Kennicutt Law (Schmidt 1959l lKennicuttl 
1998a). Feedback from supernovae is c aptured by an effective 
equatio n of state for star-forming gas ( Spri ngel & Hernquisll 
l2003al) . A prescription for supermassive black hole growth 
and feedback is also included, where black holes are rep- 
resented by collisionless "sink" particles that interact gravi- 
tationally with other components and accrete gas from their 
surroundings. The accretion rate is estimated from the local 
gas density and sound speed using a spherical Bondi-Hoyle 
feondi&HovTelll944t iBondill 19521) model that is limited by 
the Eddington rate. Feedback from black hole accretion is 
modeled as thermal energy injected into the surrounding gas 
(ISpringel et al.ll2005at iDi Matteo et al.ll2005l) . We note that 
implementations of our model for black hole growth and feed- 
back that do not explicitly account for Eddington-limited ac- 
cretion achieve similar results to the method employed by 
us (e.g. compare the works of IDi Matteo et al.l 120071 and 
ISiiacki et al.ll2007l) . 

These hydrodynamic simulations adopted the ACDM 
model with cosmological parameters from the first year 
Wilkinson Microw ave Anisotropy Probe data (WMAP1, 
ISpergel et alJl2003l) . (ft m , n h , Cl A , h, n s , <j 8 )= (0.3, 0.04, 0.7, 
0.7, 1, 0.9). In this paper, we use the same parameters. 

2.1.2. Hierarchical Assembly of the Quasar System 

In the simulation analyzed here, the quasar host galaxy 
builds up hierarchically through seven major mergers of gas- 
rich progenitors between z = 14.4-6.5. Gravitational inter- 
actions between the merging galaxies form tidal tails, strong 
shocks and efficient gas inflows that trigger intense starbursts. 
The highly concentrated gas fuels rapid black hole accretion. 
Figure Q] shows the evolution of some aspects of the system. 
Between z ~ 14-9, the merging galaxies are physically small 
and the interactions occur on scales of tens of kiloparsecs. 
By z ~ 9-7.5, when the last major mergers take place, the in- 
teractions have increased dramatically in strength. Galaxies 
are largely disrupted in close encounters, tidal tails of gas and 
stars extend over hundreds of kiloparsecs, and powerful bursts 
of star formation are triggered, resulting in an average star 
formation rate of ~ 10 3 Moyr -1 that peaks at z ~ 8.5. During 
this phase, the black holes are heavily obscured by circumnu- 
clear gas. The luminosity from the starbursts outshines that 
from the accreting black holes. So, we refer to this period 
(z ~ 14-7.5) as the "starburst phase." 

Once the progenitors have coalesced, the multiple SMBHs 
from the galaxies merge and grow exponentially in mass 
and feedback energy via gas accretion. During this period 
(z ~ 7.5-6, hereafter referred to as "quasar phase"), the black 
hole luminosity outshines that of the stars. At redshift z ~ 6.5, 
when the galaxies coalesce, the induced high central gas den- 
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FIG. 1. — Evolution of the star formation rate, black hole accretion rate, 
and masses of the black holes a nd stars, respect ively, of the simulated quasar 
system at z ~ 6.5, adopted from Li et al. (2007). 



sities bring the SMBH accretion and feedback to a climax. 
The black hole reaches a mass of ~ 2 x 10 9 M Q , and has 
a peak bolometric luminosity close to that of J 1148+5251. 
Black hole feedback then drives a powerful galactic wind that 
clears the obscuring material from the center of the system. 
The SMBH becomes visible briefly as an optically-luminous 
quasar similar to Jl 148+5251. Once the system relaxes, the 
SMBH and the host satisfy the relation, M RH w 0.002 M star 
similar to that measured in n earby galaxies (Ma gorrian et al.l 
U998tlMarconi & Huntll2003l). as a result of co-eval evolution 
of both components JLi et al] 120071: iRobertson et all 120061: 
Hop kins et aLll2007d) . 

After z < 6 (the "post-quasar phase"), feedback from star 
formation and the quasar quenches star formation and self- 
regulates SMBH accretion. Consequently, both star formation 
and quasar activity decay, leaving behind a remnant which 
rapidly reddens. The object will eventually evolve into a cD- 
like galaxy by the present day. (For an overview of this sce- 
nario, see e.g., iHopkins et al.1l2007allR) 

The photometric calculations by IRobertson et al.l d2007l) 
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FIG. 2. — Test run of the radiative equilibrium algorithm adopted from 
Bjorkman & Wood 12001) on 3-D Cartesian coordinates for a simple galaxy 
with a bulge and a disk. The input spectrum of the stars is a black- 
body with a temperature of 15000 K (black curve). The red curve is our 
output, the blue curve is simulation data from Kenny Wood, while the 
crosses are observations of a lenticular galaxy NGC 5866, which has an 
unusual extended dust disk seen exactly edge-on. Both Wood's data and 
the observations come with the code release package from Kenny Wood 
I http://www-star.st-and.ac.uk/~kw25/research/montecarlo/gals/gals.html I. 



show that this quasar system satisfies a variety of photometric 
selection criteria based on Lyman-break techniques. The mas- 
sive stellar spheroid descended from these z ~ 6 quasars could 
be detected at z ~ 4 by existing surveys, while the galaxy pro- 
genitors at higher redshifts will likely require future surveys 
of large portions of the sky (> 0.5%) at wavelengths A > 1 pm 
owing to their low number densities. 

2.2. ART 2 : All-wavelength Radiative Transfer with Adaptive 
Refinement Tree 

ART 2 is based on a unification of th e 3-D Monte Carlo 
radiati ve e quilibrium code develo ped by iBjorkman & Woodl 
d2001l) and lWhitnev et all d2003bl) . an adaptive grid, a multi- 
phase ISM model, and a supernova-origin dust model. Below, 
we describe our implementation of ART 2 . 

2.2.1. Monte Carlo Radiative Transfer for Dust in Radiative 
Equilibrium 

The Monte Carlo RT method follows the propagation, 
scattering, absorption, and reemission of "photon packets" 
(groups of equal-energy, monochromatic photons that travel 
in the same direction), by randomly sampling the various 
probability distribution functions that determine the optical 
depth, scattering angle, and absorption rates. For the dusty 
environments we are concerned with, the re-emitted spectrum 
depends on the temperature of the dust, which is assumed 
to be in thermal equilibrium with the radiation field. In the 
traditional scheme, the frequencies of the re-emitted photons 
are sampled from local reemission spectra with fixed temper- 
atures. The dust radiative equilibrium is ensured by perform- 
ing the Monte Carlo transfer iteratively until the dust temper- 
ature distribution converges. Such methods often require a 
large number of iterations and are therefore computationally 



e xpensive dBiorkman & Woodll200ll : iPascucci et aT]|2004l) . 

IBjorkman & Woodl (1200 ll) proposed a solution to this prob- 
lem by formulating an "immediate reemission" algorithm, in 
which the dust temperature is immediately updated upon ab- 
sorption of a photon packet, and the frequencies of re-emitted 
photons are sampled from a spectrum that takes into account 
the modified temperature. The advantage of this approach is 
that dust radiative equilibrium and the radiative transfer so- 
lutions are obtained simultaneousl y without iteration. Thi s 
algorithm is described in detail in Bjork man & Wood ((2001); 
here we briefly outline the steps. 

Assuming that each photon packet carries an energy £ 7 , and 
after absorption of Ni packets in the i—th grid cell, the total 
energy absorbed in the cell is 



: NiE~ 



(1) 



In radiative equilibrium, this energy must be re-radiated, 
with a thermal emissivity of j v = n u pB v {T), where k„ is the 
absorptive opacity, p is the dust density, and fi„(T) is the 
Planck function at temperature T, 



BJT)- 



2h P 



e h p u/(kT)_l 



(2) 



where hp is Planck's constant, c the speed of light, and k is 
Boltzmann's constant. The emitted energy in the time interval 
At is 



£, em = 47rAf / dViJ pK„B„(T)dv 

= 4i:AtKp(T i )B(T i )m i , (3) 

where Kp = J n u B v dv/B is the Planck mean opacity, B = 
aT 4 /ir is the frequency integrated Planck function, and m, 
is the dust mass in the cell. 

Equating the absorbed (fl]i and emitted © energies, we ob- 
tain the dust temperature as follows after absorbing N, pack- 
ets: 

*lf = NiL , (4) 
1 4NjK P (Tdmi 

where iV 7 is the total number of photon packets in the simu- 
lation, and L is the total source luminosity. Note that because 
the dust opacity is temperature-independent, the product 
Kp(Tj)aTj 4 increases monotonically with temperature. Con- 
sequently, 7} always increases when the cell absorbs an addi- 
tional packet. 

The added energy to be radiated owing to the temperature 
increase AT is determined by a temperature-corrected emis- 
sivity A j v in the following approximation when the tempera- 
ture increase, AT, is small: 



Aj v w n u pAT 



dBJJ) 
dT 



(5) 



The re-emitted packets, which comprise the diffuse radiation 
field, then continue to be scattered, absorbed, and re-emitted 
until they finally escape from the system. This method con- 
serves the total energy exactly, and does not require any iter- 
ation as the emergent SED, vL u = n„B JT), corresponds to 
the eq uilibrium temperature distribution (Bjork man & Woodl 
l200l . 

This Monte Carlo radiative equilibrium code works as fol- 
lows: first, the photon packets are followed to random inter- 
action locations, determined by the optical depth. Then they 
are either scattered or absorbed with a probability given by 
the albedo (a = n s a s /{n s a s + n a a a ), where n and a are number 
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density and cross section for either scattering or absorption, 
respectively). If the packet is scattered, a random scattering 
angle is obtained from the scattering phase function. If in- 
stead the packet is absorbed, it is reemitted immediately at a 
new frequency determined by the envelope temperature, us- 
ing the algorithm described above. After either scattering, or 
absorption plus reemission, the photon packet continues to a 
new interaction location. This process is repeated until all the 
packets escape the dusty environment. Upon completion of 
the Monte Carlo transfer, the code produces emergent SEDs 
and images at any given inclination for a wide range of broad- 
band filters, including those of Hubble Space Telescope (NIC- 
MOS bands), Spitzer (IRAC and MIPS bands), and SCUBA 
s ubmillimeter bands. 

iBiorkman & W ood (2001) tested this algorithm extensively 
by comparing a s et of benchmark calcula tions to those of 
llvezic etail d 1997b for spherical geometry. iBaes et all d2005l) 
critically study the f reque ncy distribution adjustment used by 
Bjorkman & Wood! (1200 ll) . and give a firm theoretical basis 
for their method, although it may fail for small dust grains. 

One drawback of this code, however, is its fixed geome- 
try and limited dynamic range. It can handle only 2-D or 3- 
D spherical -polar grids, which are not suitable for capturing 
the arbitrary geometry and large dynamic ranges characteris- 
tic of merging galaxies, which have multiple density centers, 
and where gas and stars extend hundreds of kpc and shocks 
produce highly condensed gas on scales of pes. We have de- 
veloped an improved version of this code by implementing an 
adaptive Cartesian grid on top of the code released by Wood 1 . 

To ensure that our implementation of the algorithm is cor- 
rect, we rerun the test problem that comes with the code re- 
lease package of Wood using a uniform Cartesian grid. The 
test problem consists of a simple galaxy with bulge and disk 
components. As is shown in Figured our code reproduces 
the result of Wood very well. 

2.2.2. Adaptive-mesh Refinement Grid 

The SPH simulations using GADGET2 dSpringell [2005b 
output hydrodynamic information as particle data. However, 
the ray tracing in the RT calculation is done on a grid. There- 
fore, it is necessa ry to interpolat e the particle-based density 
field onto a grid. Jonsson] d2006l) performed radiative trans- 
fer calculations on SPH simulations of galaxy mergers using 
an adaptive grid as implemented in his Monte Carlo RT code, 
SUNRISE. Unfortunately, self-consistent calculations of dust 
radiative equilibrium and emission are not yet included in 
SUNRISE. 

Our algorithm for c onstructing adaptive grids is similar to 
that of lJonssonl d2006l) . We typically start with a base grid of 
a 4 3 box covering the entire simulation volume. Each cell is 
then adaptively refined by dividing it into 2 3 sub-cells. The 
refinement is stopped if a predefined maximum refinement 
level, RL, is reached, or if the total number of particles in 
the cell becomes less than a certain threshold, whichever cri- 
terion is satisfied first. The maximum refinement level used 
in the present work is 12, and the maximum particle number 
allowed in the cell is 32, half the number of the SPH smooth- 
ing kernel neighbors used in the GADGET2 simulations. The 
resolution of the finest level is therefore L m ; n = Lbox/2 (RL+1) , 
where L\, ox is the box length, and RL is again the maximum 
refinement level. For example, for the parameters used in this 
simulation, L\, ox = 200 kpc, RL = 12, we have L m ; n = 24.4 pc. 

1 http://www-star.st-and.ac.uk/~kw25/research/montecarlo/gals/gals.html 



The adaptive-mesh refinement grid serves as an efficient 
tree for mass assignment owing to fast neighbor finding within 
the grids for SPH smoothing. After the grid is constructed, 
the gas properties at the center of each grid cell, such as 
density, temperature, an d metallicity, are calcula ted using the 
SPH smoothing kernel dHernquist & Katall989h of the origi- 
nal simulation. All physical quantities are assumed to be uni- 
form across a single cell. 

Figure [3] gives an example of the adaptive grid a pplied to 
a snap shot from the galaxy merger simulations in iLi et alj 
(2007), and compared with a uniform grid. This particu- 
lar snapshot represents the time when the system reaches the 
peak quasar phase, when the galaxies are in the final stages of 
coalescence. The system is highly dynamical as much gas is 
falling into the center, while feedback from the central mas- 
sive black hole drives an outflow. The gas distribution is thus 
inhomogeneous. 

As is apparent in Figure[3] a uniform grid of 50 3 (top panel) 
barely captures the density distribution with a spatial resolu- 
tion of 4 kpc. The resolution is worse than an adaptive grid 
with a moderate maximum refinement level (i.e., RL=5, mid- 
dle panel), which has a finest cell length of ~ 3.1 kpc. The 
grid with a maximum refinement level of 12 shown in the bot- 
tom panel fully captures the large dynamic range of the gas 
density distribution in three dimensions. It has a minimum 
cell length of ~ 24.4 pc for the finest cells, which is compa- 
rable to the spatial resolutio n of the original hydrodynamic 
simulation of ILi et all d2007l) . This resolution is equivalent to 
a ~ 10000 3 uniform grid, which is impractical with existing 
computational facilities. 

2.2.3. A Multiphase Model for the Interstellar Medium 

In determinin g the dust distribution , we ad opt the multi- 
phase model of Springel & Hernquist (2003a) for the ISM. 
The ISM is then comprised of condensed clouds in pres- 
su re equilibrium with an ambient hot gas, as in the picture 
of iMcKee & Ostrikerl d!977l) . In the hydrodynamic simula- 
tions, individual SPH particles represent regions of gas that 
contain cold, dense cores embedded in a hot, diffuse medium. 
The hot and cold phases of the ISM co-exist in pressure equi- 
librium but have different mass fractions and volume filling 
factors (i.e., the hot-phase gas is < 10 % in mass but > 99% 
in volume). In our previous studies, ILi et alJ (2007) and 
iHopkins et al. (2006) used only the hot-phase density to de- 
termine e.g. the obscuration of the central AGN, as the major- 
ity of sight lines will pass through only this component owing 
to its large volume filling factor. However, this method gives 
only an effective lower limit on the column density. 

Here, we consider two components of the dust distribution, 
having different dust-to-gas ratios and being associated with 
the two phases of the ISM. Within each grid cell in the RT 
calculation, the hot gas is uniformly distributed, while the 
cold, dense cores are randomly embedded. Because of the 
much higher density and smaller volume of the cold phase, it 
is impractical to resolve these clouds either in hydrodynamic 
simulations or in our radiative transfer calculations. Con- 
sequently, we implement an observationally-motivated, sub- 
resolution prescription to treat the cold clouds, constrained by 
the observed mass spectrum and size distribution of molecular 
clouds in galaxies. 

Observations of giant molecular clouds (GMCs) in g alaxies 
show that the GMCs follow simple scaling relations dLarsonl 
1 198 lb . namely a power-law mass distribution, d/Y/dM oc 
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FIG. 3. — Example of the adaptive grid applied to a snapshot of the galaxy merger from Li et al. (2007). From top to bottom is: uniform grid, and adaptive 
grids with maximum refinement levels of RL=5, and RL=12, respectively. The label GN indicates the total grid number covering the entire volume (in case of 
adaptive grids, it provides the number of the finest grid and that of the total grid). This plot demonstrates that adaptive grids with sufficient refinement capture 
the density distribution well, while the uniform grid fails to do so unless an unreasonably large number of cells is employed . Throughout this paper, we use a 
maximum refinement level of 12, which has a grid resolution comparable to that of the hydrodynamic simulations in Li et al. 12007). 
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M 2 (e.gjFuller & Mversl 19921:1 Ward-Thompson et al.lll994t Integrating over the cloud radius, one obtains 



lAndre et alJ 119961: iBlitz & Rosolowskyl 120061) , as well 



as a 



powe r -law mass-radius r elation M oc R~ (.e. g, Sanders et al. 
1985 1: iDame et alj|1986fc IScoville et afll 19871 ISolomon et al.l 
1987t lRosolowsky l l2005L 12007) . It has been suggested by the- 
oretical modeling that the mass funct ion is produced by turbu- 
lence i n self-gravita t ing cl o uds (e.g..lElm egreen & Falgarone 
1996 1 : lElmegreenl 20021 : [Ballesteros-Paredes & M ac Low 
20021: iLi et al] 120031 while the mass-radius relation is at- 
tributed to virial equilibrium (Larson 1981b. 

Here we incorporate these two empirical relations into our 
ISM model for the RT c alculations. For a given cell , the 
two-phase break-down in Springel & Hernquist (20033) de- 
termines the hot and cold phase gas density according to pres- 
sure equilibrium. The cold clouds are assumed to follow the 
Larson scaling relations: 



dn 



(6) 



■ =AM~ 
dM 

M = BR , (7) 
where dn / dM is the number density of the clouds differen- 
tial in cloud mass M, and R is the cloud radius. From these 
equations, one obtains the cloud size distribution 

— =/3AB 1 - a R< a ^ 1 -^ 



dR 



(8) 



where C = [3AB l ~ a , and 7 = a(3+ 1 - (3. 

Assume that the minimum and maximum values of the 
cloud mass are Mq and M\ , then the normalization constant 
of the mass spectrum for each cell can be determined using 
dn 

--Xcpc, (9) 



or 



dM 



■■X c p t 



-MdM ■ 



M\~ 



-M\; a ' 



(10) 



where p c and x c are the cold gas density and volume filling 
factor, respectively. The normalization constant of the M-R 
relation may be determined using 



4ir o dn 



-R" 



dM 



dM, 



or 



where 



4-nA 



3r)x c 



0/3 



1 + 



(ii) 



(12) 



(13) 



The minimum and maximum cloud radii in the cell are there- 
fore 



^0= ~TT 



Mo 
B 



Ri 



V/3 



V/3 



(14) 



For a photon traveling a distance L in the cell, the average 
number of cold clouds of radius R the photon will intersect is 
given by 



dN n r dn 
— = ttR L — 
dR dR 

= nCR 2 -~ f . 



(15) 



R 



N = itLC- 



3-7 p3-7 

3-7 



(16) 



Therefore, the average distance the photon must travel to hit a 
cold cloud (the mean free path) is given by 



"7 



ttC 



(r^-r^) 



(17) 



In the dust RT calculation, we assume that dust is associated 
with both the cold and hot phase gases through certain dust- 
to-gas ratios. When a photon enters a cell, we first determine 
the distance L/, it travels in the hot phase gas before hitting a 
cold cloud, which is an exponential distribution function 



p(L) = — exp 



L 

Li,, 



(18) 



Therefore, L/, = -L„,ln£, where £ is a random number uni- 
formly distributed between and 1 . The radius of the cloud 
the photon just hits is also determined randomly assuming the 
distribution function of Equation ( fT3l >: i.e., 



R = 



p-i-7 
^0 



1/(3-7) 



(19) 



The distance L c the photon travels in this cold cloud is again 
a random variable given by L c = 2R^, because clouds are 
assumed to be uniformly distributed. These equations com- 
pletely define the statistical procedure for determining the 
dust column densities associated with Lf, and L c as: 



■■ PhUi 
3BR> } - 



4tt 



(20) 



With a given dust opacity curve, these equations allow one 
to calculate the optical depths, 17, and t c for the hot and cold 
dust, respectively. They relate the photon path lengths in the 
multiphase ISM to the total optical depth r tot , which is then 
compared with a randomly drawn number 77, to determine 
whether the photon should be stopped for scattering or ab- 
sorption. In detail, the Monte Carlo ray-tracing procedure for 
the radiative transfer therefore involves the following steps: 

1 . A photon packet is emitted from either a blackhole or a 
stellar source with random frequencies consistent with 
the source spectra. The photon is emitted with a uni- 
formly distributed random direction. The probability of 
a photon being emitted by any given source is deter- 
mined by its luminosity relative to the total. 

2. A random optical depth over which the photon must 
travel before an interaction with the dust occurs, t, = 
-ln£, is drawn to determine the interaction location. 
The interaction includes scattering and absorption. In 
our method, the photon energies are not weighted, only 
one event is allowed. That is, at any given interaction 
site, the photon is either scattered or absorbed, but not 
both. 

3. Starting from the location of the photon emission, the 
cumulative optical depth of the photon, T tot , is calcu- 
lated stochastically using Equationl20lfor both hot and 
cold dusts. If the photon is stopped for interaction 
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within a single cell, then r to t is the sum of contribu- 
tions from possibly multiple segments of both hot and 
cold dusts within this cell. If the photon passes through 
multiple cells before an interaction occurs, then r tot is 
the sum of all contributions from relevant segments in 
these cells. 

4. At each boundary between the hot and cold phase gas 
clouds, or at the boundary of the grid cell, the next 
interaction point is determined by the comparison be- 
tween t, and T tot . If t, < T tot , then the photon is ei- 
ther scattered or absorbed, with a probability given by 
the scattering albedo. The exact interaction location 
is then determined inside either hot or cold phase gas, 
such that rtot becomes exactly 77. If the photon is scat- 
tered, its direction and polarization state are altered 
using the Henyey-Greenstein phase function, and the 
ray-tracing of the new photon is repeated from step 2. 
If the photon is absorbed, depending on whether the 
absorption site occurs in the hot or cold phase dust, 
the temperature of the appropriate dust cell is raised, 
an d a new photon is re e mitted according to the scheme 
of iBjorkman & Woo d (2001). The ray-tracing of the 
newly emitted photon again restarts from step 2. 

5. If the photon escapes from the system without reaching 
the optical depth r,, it is then collected in the output 
spectrum and image. The next photon will be picked up 
from the source, and the whole Monte Carlo procedure 
from step 1 will be restarted. 

After all N 1 photons have been traced, one obtains the dust 
temperature distribution for both hot and cold dusts in radia- 
tive equilibrium, as well as the output spectra and images. 
Note that such a stochastic procedure outlined here may not 
be as accurate as physically tracking the locations and sizes 
of the cold phase gas clouds, which can be rather difficult, if 
not impossible. However, this method works efficiently for a 
large number of cold clouds uniformly distributed in the cells, 
and it gives correct average extinction properties for a multi- 
phase ISM model we are interested here. 

Hereafter, we refer to "HPG-dust" as the dust that originates 
from the hot-phase gas, and "CPG-dust" as that originating 
from the cold-phase gas. Note that the gas only determines 
the distribution and mass of the dust through a given dust-to- 
gas ratio. The dust temperature is not associated with the gas 
temperature, and is calculated self-consistently according to 
radiative equilibrium. We further refer to "cold dust", "warm 
dust" and "hot dust" as dust with temperatures T < 100 K, 
100 < T < 1000 K, and 1000 < T < 1200 K, respectively. 

In the RT calculations, we adopt a mas s spectrum with a = 
1.8, as suggested by the observations of Blit z & Rosolowskvl 
(2006), and an observed mass-radius relation with = 2.0, 
which is also a result of the virial theorem. The resulting 
mass-radius relation in our simulations has a normalization in 
the range ~ 10- 10 4 M Q /pc 2 . For a cloud size of ~ lpc, the 
normalizatio n is ~ 300MfT)/pc 2 , similar to observatio ns of the 
Milky Way dScoville et alj|1987t lRosolowskvll2007t) . It has 
been shown that the normalization of the mas s-radius relation- 
ship depends on g alactic environment (e.g.,[ElmegreeiJT989|; 
lRosolowskvll2005l:lBlitz et alj|2007l) . For example, it is about 
~ 50 M Q /pc 2 in the Large Magellanic Cloud but could be two 
orders of magnitude higher in ULIRGs. In extreme starburst 
galaxies, the normalization could go up to 10 4 M Q /pc 2 , as we 
find in our simulations. 
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FIG. 4. — Comparison of the dust absorption opacity curves from a 
supernova-origin dust model (SN model, with silicate fraction f s = 10%, 20% 
in red and black, respectively) and those of Weinaartner & Draine 12001) 
(WD01 model, with Ry = 3.1, 4.0, 5.5, in blue, cyan, and green, respec- 
tively). Note there are two differences between these two models: the silicate 
feature at ~ 9.7 /im, and a higher opacity in the optical band. 



We assume the cold clouds to be in the mass range of 
M = 1O 3 M and M\ = 1O 7 M , which is similar to that of 
protoclu s ters c l ouds in s tar fo rming galaxies, as shown in 
iLietalJ d2004 l2005allbl 120061) . For the cold clouds, we 
have enhanced the pressure by a factor of 10 over the ther- 
mal pressure to account for the effects of turbule nce (e.g., 
iBlitz & Rosolowskvl2006tlChakrabarti et alj2007al) . This en- 
hancement factor is referred to as "CP" hereafter. The dust- 
to-gas ratio of the cold clo uds is chosen to be the same as 
the Milky Way value (l:124; IWe"ingartner & Draindl2001l). as 
found in a large sample of ULIRGs (Dunne & Eales 120011 : 
iKlaas et ai1l2001h . while that of the hot, diffuse gas is chosen 
to be 1 % of the Milky Way value, consistent with the dust 
survival rate of sputte ring in a hot, diffuse ISM dBurke & Silk! 
1 1974t [Reynolds et alj|1997l) . Observations of some obscurred 
or red AGNs also su ggest a ratio significantly lower than that 
of Milky Way (e .g., Maiol mo et alj200"ltlKuraszkiewicz et alj 
l200llHallet"al . 2006). We will perform a systematic param- 
eter study of these choices in §|4] However, we note that these 
values are found to best reproduce the observed quasar SED 
at z ~ 6.5. Because the dust opacity is proportional to the 
gas metallicity, therefore the dust opacity is weighted by the 
metallicity of the gas for both hot and cold gas phases. 

2.2.4. Supernova-origin Dust Model 

Our understanding of the formation and distribution of dust 
has benefited from dust m aps of the Milk y Way (e.g.. | Gehrzl 
119891: ISchlegel et al.|[l998h . In particular, |Gehrz| dl989j) con- 
cludes from a detailed Galactic survey of dust-producing stars 
that dust in the Milky Way originates in three principle ways: 
(1) by condensation in winds of evolved, post-main-sequence 
objects, which accounts of ~ 90% of the stellar dust; (2) by 
condensation in ejecta from massive novae, supernovae and 
Wolf-Rayet stars, which amounts to < 10%; and (3) b y slow 
accretion in molecular clouds (see also Marchenko 200g for a 
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FIG. 5. — Distribution of the optical depths of both HPG- and CPG-dust, 
respectively, at A = 0. 1 fim where the dust opacity peaks as shown in Figure[4] 
The density grid corresponds to the simulation snapshot at z = 6.5 shown in 
Figure [3](bottom panel). The "HPG-dust" and "CPG-dust" refer to the dust 
a ssociat ed with hot-phase and cold-phase gas, respectively, as described in 

sum 



review). 

Over the past several decades, various dust models have 
been developed based on the observed extinction curves 
of the Large Magellanic Cloud (LMC), the Small Magel- 
la nic Cloud (SMC) and the Milky Way (e.g., see reviews 
bvlSavage & Mamis||197^lMathislll99"otlCalzetti et al.lll994t 
Dorschner & Hennind 119951: ICalzetti et al.1 120001: IWhittetl 
20031: iDraind 120031) . In these classical models, dust is 
assumed to form in the envelopes of old, low-mass stars 



such a s asymptotic giant branch (AGB) stars with ages > 1 
Gvrs (lMathis|l l990: Mor gan & Edmundsll200l lDwell200l 
lMarchenkol2006l) . 

However, this picture may be different for young, high- 
redshift objects. Recent deep millimeter and submillime- 
ter observations of several z ~ 6 quasars, which trace the 
far-infrared thermal dust emission from these systems, show 
large masses of dust in these quasar hosts when the Uni- 
verse was less than 1 Gyr old (e.g., [B ertoldi et al. 2003a; 
Char mandaris et al.1 120041 iRobson et al.l l2004t ICarilli et al l 
2004t i Maiolino et al J 12004 iBeelen et all 120061: I ffines et al.1 
2006t Ijiang et al.l l2006t IWillott et al.l 120071) . In particular, 
Maiohno et al.l (120041) find that the extinction curve of the 
reddened quasar SDSS J1048+46 at z = 6.2 is different from 
those obser ved at z < 4 (similar to that of the Small Magel- 
lanic Cloud, Hopkins et al. 2004), but matches the extinction 
curve expected for dust produced by supernovae. 

It has been suggested that core-collapse supernovae (Type- 
II SNe) may provide a fast and efficient mechanism for 
dust production. The observational evidence for dust for- 
mation in SNe comes from observations of SN1987A (e.g. , 
Gehrz & NevlfT987t iMoselev et al.lfl989l: iRoche et al.l ll993h 
Spvromilioet al.1 119931: IColgan et all 1 19941) , Cassiopeia A 
(iDunne et al.l 120031: iDwekl l2004t see however iKrause et all 
2004 who argue that the dust emission is not associated 
with the remnant), Kepler's supernova remnant (Morg arTet al.l 



2003), and SNe 2003gd dSugerman et al.l 12006). Theoreti 
cally, several groups have calculated dust fo r mation in the 
ejecta of Type-II SNe (fTodini & Ferraral 12001 1: iNozawa et al ' 
| 2003[:ISchneider et al.l2004tlHirashita et al .1120051: iDwek et al 
120071) . In particular. iTodini & Ferraral (1200 ll) have developed 
a dust model based on standard nucleation theory and tested 
it on the well-studied case of SN1987A. They find that SNe 
with masses in the range of 12-35M Q produce about 1% of 
the mass in dust per supernova for primordial metallicity, and 
^3% for solar metallicity. 

In the present work, we adopt the dust size distribution of 
ITodini & Ferraral (1200 ll) for solar metallicity and a M = 22 M Q 
SN model, as in Figure 5 in their paper. This size distribution 
is then com bined with the dust abso r ption and scattering cross 
sections of IWeingartner & Draind d2001l). to calculate dust 
absor ption opacity curves dFinkbeinerlll 999; Finkbein er et alJ 
1999). We note that Bianchi & Schneider (2007) re-analysize 
the model of ITodini & Ferrara nfeOOll) and follow the evolu- 
tion of newly condensed grains from the time of formation to 
their survival. This new feature shows better agreement with 
observations and further supports our motivation to use a SN 
dust model. Figure [4] shows dust absorption opacity curves, 
in the range 10~ 2 - 10 4 /im, from bot h the supernova-origin 
dust m odel (hereafter SN model) and Weingartn er & Draind 
d2001l) (hereafter WD model), respectively. The SN models 
include silicate fractions / s = 10% (in red) and / s = 20% (in 
black), which show slight differences in the silicate feature at 
~ 9.7 fim (/ s = 20% model has a slightly higher opacity). The 
WD models include all the curves commonly used with ex- 
tinction Ry = A ( B yl(y) =3.1, 4.0, 5.5, which differ in the opac- 
ity in the UV-NIR (0.01 fim- 1 ^m) bands (R y = 3.l model has 
the highest opacity). 

Between the SN and WD models, there are two noticeable 
differences; one is the silicate feature at ~ 9.7 /im, and the 
second is the UV-NIR band. On one hand, the WD model has 
strong peaks around ^9.7 fim (silicate feature), while the SN 
model produces an opacity lower by nearly one order of mag- 
nitude, because the ejecta of SNe with mass above 15 M© pre- 
dominantly form amorphous carbon grains, which decreases 
the silicate fraction in the dust grains. On the other hand, the 
SN model increases the opacity in the optical band by a fac- 
tor of a few owing to the smaller grain size. We note that the 
small grain sizes may cause quantum fluctuations in the tem- 
peratu re and high-temperature transients (e.g., iDraine & Lil 
120011) . However, the dust produced by the SN model is dom- 
inated by graphite grains, which have a size distribution from 
~ 100-2000 A, comparable to the size range of the dust 
grains in the WD model, so the temperature fluctuation is ex- 
pected to be insignificant. Note we do not include polycyclic 
aromatic hydrocarbon (PAH) features at ~ 7.7 u rn, which was 
model ed in detail by iLi & Draind (120021) and IDraine & Lil 
J2007I) . The PA H feature is a goo d diagnostic of starbursts 
at low redshifts dHouck et al.l l2005). however it is very diffi- 
cult to detect in z > 6 systems. Also dust spin is not included 
in our modeling, which could be significant aro und 10 mm 
dDraine & Lazarianlll998l:lFinkbeiner et al.ll2002l) . 

The differences between these two models, in particular 
the silicate feature, may be diminished by c hoosing differ- 
ent gr ain compositions and size distributions (lLaor & Draind 
1993). For example, the WD model would converge to the 
SN mod el if the silicate fraction is reduced to ~ 5%. We 
note that lElvis et al.l (120021) propose an alternative mechanism 
for dust production in quasars. They suggest that the physi- 
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cal conditions (i.e., low temperature and high density) in the 
quasar outflow are similar to those in the envelopes of AGB 
stars, and hence may produce dust similar to that made by 
AGB stars. This scenario does not have the same timescale 
issue as does the AGB model, but only requires heavy metals 
such as carbon and silicate, which should be available even 
in high-redshift quasars, as noted in the introduction. There- 
fore, quasar winds may also serve as efficient factories for 
producing dust with extinction curves similar to those of the 
WD model. Currently there are no observations available to 
distinguish between these models at high redshifts. We will 
return to this in § 14.41 

Figure [5] shows an example of the resulting optical depths 
of both HPG- and CPG-dust at wavelength A = 0. 1 /im where 
the dust opacity reaches its maximum, as shown in Figure [4] 
This plot demonstrates that the optical depth is resolved very 
well in the hot, diffuse gas which occupies ~ 99% of the 
volume, while the optical depth of the cold, dense cores is 
usually much larger than unity, in particular in optical bands 
with high frequencies. The energy absorbed by the cold dust 
in these wavelengths is then re-emitted at infrared or longer 
wavelengths. 

3. THE SPECTRAL ENERGY DISTRIBUTION OF A QUASAR 
SYSTEM AT Z ~ 6.5 

We now calculate the spectral energy distributio n from 
UV/op tical to submillimeter of the modeled quasar in lLi et af] 
d2007l) by applying ART 2 to our SPH simulations, which pro- 
vide the gas density field and heating sources for the RT calcu- 
lations. The input spectrum includes that from stars and black 
holes, as shown in Figure [6] (top panel). The stellar spec- 
trum is calculated using the stellar population synthesis code 
STARBURST99 (iLeitherer et al]fl999t IVazauez & Leithereil 
2005). The age, mass and metallicity of the stars a re taken 
from the SPH simulations of quasar formation in iLi et alj 
d2007l) . while the stellar initial mass function (IMF ) is as- 
sumed to be a top-heavy Kroupa IMF (lKroupall2002l) . which 
characterizes a starburst better than the classical Salpeter IMF 
(ISafoeterill955l) . 

The input black hole spectrum is a composite template from 
iHopkins et alJ (l2007d|). whi ch consists of a broken power-law 
(e.g.. lLaor & Draindll99l iMarconi et all 120041) and an IR 
component thought to come from the hottest dust around the 
AGN, which cannot be fully resolved in our hydrodynamic 
simulations. The normalization of this spectrum is the to- 
tal bolometric luminosity of the black hole s. Compared to 
the luminosity calculation in ILi et al.l (|2007), the black hole 
spectrum here is normalized to the bolometric luminosity 
of Jl 148+5251. Heating by cosmic microwave background 
radiation at different redshifts z is also taken into account 
by including a uniform radiation field with temperature of 
r cmb = 2.73x(l+z)K. 

The emergent spectrum has a wavelength range of 10 2 - 
2 x 10 7 in the rest frame, and 50 viewing angles (10 in polar 
angle evenly divided in cos6> and 5 in azimuthal <f>). We use 
10 7 photon packets isotropically emitted from sources, and 
the maximum refinement level of the adaptive grid is RL=12, 
which is above the requirement for convergence (see § 14. ll for 
resolution studies). 

The calculated rest-frame SEDs of the system during the 
peak quasar phase at z = 6.5 are shown in Figure [6] (bottom 
panel). The output SED depends on the viewing angle. The 
red curve is the isotropically averaged SED (e.g., it can be 
understood as an average SED multiplied by the solid an- 



gle 471"). This average d SED agrees ver y well with observa- 
tions of Jl 148+5251 dJiang et al.1 12006). A prominent fea- 
ture of this SED is the two infrared bumps peaking around 
3 /im and 50 /im. The 50 /im bump is produced by cold dust 
(T ~ 50 K) heated by strong star form ation, as commonly 
seen i n the SEDs of starburst galaxi es (ISanders & M irabel 
Il996t ISiebenmorgen & Kriigell 120071) . The 3 /tm bump is 
produced by the hot dust (T ~ 1000 K) heated by the cen- 
tral AGN. This unique feature appears to be ubiquitous in 
most of the quasar SEDs over a wide range of redshifts (e.g. , 
l Elvis et al.lll994t IVanden Berk et al.ll200U iTelfer et alj|2002t 
Vigna li et alJl2003l:lRichards et alJl2006HJiang et alJ l2006). 

To demonstrate the line-of-sight dependence, Figure [6] also 
shows the range of SEDs viewed from two orthogonal an- 
gles of the Cartesian grid: along the z-axis (upper range) 
and the xy-plane (lower range), respectively. This range in- 
dicates difference in column density along the sight lines, and 
shows that dust extinction in the UV/optical bands differs by 
a factor of ~ 3, but no difference at wavelengths longward 
of 1 /<m. This suggests that the dust is close to optically thin 
along these two viewing angles during this time, and that in- 
frared or submillimeter observations of luminous quasars may 
not be diagnostics for the orientation of the host. We have 
also checked the three major components of the output SED, 
namely the scatter, escape and reemission, and found that the 
photon energy is conserved. This confirms the photon con- 
servation algorithm used in the radiative transfer calculation. 
Moreover, we find that the emergent SED in the UV/optical 
bands (0.01 - 1 /im) is dominated by scattering. 

It is interesting to comment on the spectral feature at wave- 
length A r est = 9.7/xm where the silicate cross section peaks 
(Draine 2003). This feature can produce either emission or 
absorption in the spectrum depen ding on the optical dept h 
of the medium at this wavelength (Bio rkman & Woodll2001|) . 
In order to have absorption at this wavelength, the medium 
has to be optically thick (i.e., T^ri^m >> !)■ Generically, 
this would also result in deep absorption in the optical/NIR 
bands (~ 0.01 - 1 /im) owing to their much higher dust opac- 
ity, which is almost two orders of magnitude higher than that 
at 9.7/im according to the dust opacity curve in FigureH] Our 
spectra exhibit absorption when the system is in the starburst 
phase (e.g., z > 10) when the object is highly obscured. As the 
system proceeds to the quasar phase (e.g., z < 8), the emis- 
sion feature becomes increasingly prominent as dust becomes 
more and more transparent to the radiation. 

In observations, both absorption and emission features at 
Arest = 9.7/im hav e been detected for a wide range of ob - 
jects up to z ~ 3 (lArmus et al.ll2004t iPapovich et alJ 12006). 
An absorption feature is fre quently seen in du sty star-forming 
galaxies such as M 82 fe.g.. lSturm et alJl2000l) . Arp 220 (e.g., 
ISpoon et alJl2004l) . and ULIRG IRA S08572+39 1 5 whic h ex- 
hibits the most extreme absorption dSpoon et alj I2006T) . A 
comprehensive collecti on of the SEDs of starbur s t nucl ei and 
ULIRGs is reviewed by lSiebenmorgen & Kriigell d2007l) . who 
also provide a library of 7000 SEDs for dusty galaxies. 

On the other hand, an emission fea t ure is also reported 
in many observations dHao et al.l 12005: Siebenmorge n et all 
120051: ISturm et alj|2005l) . In particular. lHao et alJ d"2007l) an- 
alyze a sample of 196 local AGNs and ULIRGs obser ved 
by the Infrared Spectrograph (IRS; i Houck et al.l 12004) on 
board the Spitzer Space Telescope d Werner et alJ 120041) to 
study the distribution of strengths of the 9.7/im silicate fea- 
ture. These authors find a wide range of silicate strengths: 
quasars are characterized by silicate emission and Seyfert Is 
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Fig. 6. — SEDs of a quasar system at z — 6.5 determined by applying ART 2 to the SPH simulations of iLi et all d2007h. The inpu t spectra 
(top panel) include the stellar spectrum (in bluel calculated with STARBURST99 (Leitherer et al. 1999; Vazquez & Leitherer 2005), and the 
composite spectrum of the AGN (in black, Hop kins et alj |2007d). Note this snapshot corresponds to the peak quasar phase of the system, 
the stellar radiation is insignificant compared to that from the AGN at this evolution stage. The red curve is the total input spectrum. The 
output spectra are shown in the bottom panel. The red curve is the total SED assuming an isotropic distribution of photon energies in all 
directions, while the grey region indicates the range from two orthogonal an gles corre sponding to the z-axis and xy-plane of the Cartesian grid, 
respectively. The filled black circles with error bars are observations from lJiang et al.l (|2006). 
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equally by emission or weak absorption. Seyfert 2s are dom- 
inated by weak silicate absorption, and ULIRGs are charac- 
terized by strong sili cate absorption (me an apparent optical 
depth of about 1.5). I Spoon etaT] d2007t) find that the same 
sample of galaxies is systematically distributed along two dis- 
tinct branches: one with AGN-dominated spectra and one 
with deeply obscured nuclei and starburst-dominated spec- 
tra. These authors suggest that the separation may reflect a 
fundamental difference between the dust geometries in the 
sources: clumpy for AGNs versus non-clumpy o bscuration 
for starbursts. For example, iLevenson et alj d2007l) suggest 
that the extremely deep absorption in IRAS08572+3915 re- 
quires a source to be embedded in a smooth distribution of 
material that is both geometrically and optically thick. Our 
simulations show that the dust distribution around the quasar 
is quite clumpy, and that the SEDs changes from starburst to 
AGN-dominated as the system evolves from the starburst to 
quasar phases (see Figure [Toll. These are consistent with the 
observations. 

4. PARAMETER STUDIES 

In this Section, we explore the large parameter space in- 
volved in this radiative transfer calculation, and systemati- 
cally study the resolution convergence, parameters in the two- 
ISM model, input spectra for the BHs, the dust models, and 
the dust-to-gas ratios for both the cold- and hot-phase gas. 

4.1. Resolution Studies 

Figure [7] shows resolution studies for photon number (top 
panel) and grid size (bottom panel). The SEDs converge when 
the photon number is larger than 10 6 . If the photon number 
is low, then Poisson noise is significant, which results in large 
fluctuations in the SEDs. Therefore, throughout the paper, we 
use a photon number of 10 7 for SED production, and 10 8 for 
images in order to have higher signal-to-noise. 

For the adaptive grids, as the refinement level increases, 
the hot dust in the central region around the AGN is bet- 
ter resolved, contributing to the hot dust bump in 1 - 10 /im. 
The SEDs converge when the refinement level goes above 10, 
which has a minimum cell size close to the spatial resolution 
(~30 pc) of the original hydrodynamic simulations. Com- 
pared to the adaptive grid method, SEDs using uniform grids 
with reasonable computing expense have poor resolution. As 
demonstrated in Figure[5J a 50 3 uniform grid has a resolution 
of 4 kpc, which can only resolve the dust in the diffuse gas in 
the outskirts of the system, but not the clumpy, dense regions 
around the central AGN. Consequently, the resulting SEDs do 
not entirely resolve dust emission in 1 - 10 fim which comes 
partly from hot dust near the AGN, and the cold dust bumps 
are lower by up to one order of magnitude. 

4.2. ISM Model Parameters 

Figure [8] shows the SEDs with various parameters for the 
two-phase breakdown of the ISM. When there is no phase 
break-down (purple curve), only hot-phase gas is considered 
(no cold gas). In this case, the hot gas has both large volume 
filling factor and mass fraction, so the extinction is extremely 
high, which leads to significant emission in the NIR but little 
emission in the 20- 1000 /^m range (the cold dust is sparse). 
In cases with cold gas only (no hot-phase gas), the SEDs show 
big cold dust bumps but no hot dust emission, as indicated by 
the blue and cyan curves. A comparison of these two curves 
shows that a larger cold phase pressure leads to a higher vol- 
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FIG. 7. — Resolution studies of the number of photons (top panel), and 
the grid refinement level (bottom panel), respectively. The top panel shows 
the SEDs produced with photon number N„h = 10 4 — 10 9 . These SEDs have 
similar shapes, and they converge when N p h > 10 6 . A smaller photon number 
results in larger fluctuations in the SED owing to greater Poisson error. The 
bottom panel compares the SEDs produced with a uniform grid and those 
with adaptive grids of different refinement level. Compared to an adaptive 
grid method, SEDs using uniform grids with grid number GN = 30 3 and 50 3 
do not have sufficient dynamic range to resolve both the cold and hot dust, 
which result in an underestimate of the dust emission longward of 10 /im. As 
the grid refinement level increases, hot dust is better resolved, contributing 
to the hot dust bump at 1 - 10 fj,m. The SEDs converge when the refinement 
level goes above 10, which has a minimum cell size approaching that of the 
spatial resolution of the original hydrodynamic simulations. Level 12 is the 
standard refinement level used in this paper. 



ume filling factor for the cold dust, resulting in stronger cold 
dust emission. 

In the cases where both hot and cold phases coexist, the 
amount of cold dust emission depends on the cold gas pres- 
sure, as indicated by the other colored curves. We find that 
cold pressure in the range of 10 - 100 is able to produce 
cold dust bumps that fit the submillimeter observations of 
Jl 148+5251. However, because the SED with CP=10 has 
both the hot and cold dust bumps that agree better with the 
mean SEDs o f lumi nous quasars in the Sloan samples of 
iRichards et ail d2006l) . we choose CP=10 as a standard value 
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FIG. 8. — Parameter study of the two-phase break down of the ISM with hot 
and cold phases. In the legend, HP=1 indicates existence of hot-phase gas, 
while C P=10 indicates that the cold pressure is enhanced by a factor of 10 (see 
jj |2.2.3l for more details). The purple curve represents the "no phase break- 
down" case in which the cold gas is not considered; there is only hot-phase 
gas whose density is the same as that given directly by the SPH simulations. 
The blue and cyan curves represent cases in which no hot-phase gas is present 
(only cold-phase gas is considered), but the pressure enhancement factor for 
cold-phase gas varies from 10 to 100, respectively. The rest of the colored 
curves represent cases in which both hot and cold phases co-exist, with cold 
gas pressure varying from 1 to 100. In the RT calculations, we use standard 
values CP=10andHP=l. 
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FIG. 10. — A parameter study of the input spectrum for the black hole. The 
solid curves represent the output SEDs, while the dotted lines represent the 
input black hole spectra. Th e blue and red curves indicate the use of an intrin- 
sic, broken power-law as i n Marconi et al. (2004), and a composite spectrum 
from Hopkins et al. 1 20()7d), respectively. Both input spectra have the same 
bolometric luminosity. A power-law input spectrum would require pc-scale 
resolution to resolve the dust near the AGN in order to produce the near-IR 
emission, which is below the resolution of our hydrodynamic simulations. A 
composite spectrum therefore serves as a sub-resolution recipe to resolve the 
dust emission within parsecs of the AGN. 
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FIG. 9. — Parameter study of the power-law indices o f the m ass spectrum 
a and mass-radius relation /3 of the cold clouds (see § 12.2.31 for more de- 
tails). The output SED is not sensitive to the ranges of a and /3 considered 
here. In the RT simulations, we adopt a = 1.8 and /3 = 2.0 as suggested by 
observations. 



in our calculations. 

Figure [9] shows a study of the output SEDs obtained by 
varying the power-law indices of both the mass spectrum 
and the mass-radius relation of the cold clouds. Increasing 
a would steepen the cloud mass function, leading to more 



smaller cold clouds that boost the cold dust bump. A simi- 
lar effect is seen by increasing (3. However, the SEDs are not 
very sensitive to the change of either a or (3 in the ranges of 
1.5 - 2.5 and 1.5 - 3.0, respectively. Changing a from 1.5 
to 2.5, or (3 from 1.5 to 3.0 only results in a change in the 
SED by a factor of about 2. In the RT calculations, we adopt 
the standard values of a = 1.8 and (3 = 2.0, as suggested b y 
observations (Blitz & Rosolow skv 20061: lRosolowskvll2007h . 

4.3. Input Spectrum 

Figure [10] shows a comparison of the emergent SEDs using 
different input black ho l e spec trum, namely a broken power- 
law as in lMarconi et al. (2004) (blue c urve), and a composite 
spectrum from Hopkins et al. (2007d) (red curve). All other 
parameters being equal, the SED using the power-law spec- 
trum shows more emission in the optical bands and lower 
emission in the IR than the SED with the composite input 
spectrum. These differences owe their origin to differences in 
the input spectrum. The composite spectrum includes emis- 
sion from hot dust residing in the vicinity of the AGN that is 
below the resolution limit of our hydrodynamic simulations. 
This figure emphasizes the care that must be taken when in- 
cluding an AGN and performing radiative transfer on scales 
that are not well-resolved, and demonstra t es that a composite 
black hole spectrum as in iHopkins et alj (|2007d) provides a 
viable sub-resolution prescription for the dust emission con- 
tributed within parsecs of the AGN. 

4.4. Dust Models 

Figure [TTJ shows the emerg ent SEDs using dust extinc- 
tion curves from the WD model ( Weing artner & Draind200ll) 
with Ry = 3.1 and our SN model, respectively. The SEDs 
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FIG. 11. — Comparison of the SEDs from using the dust extinction curves 
of Weingartner & Draine (2001) with Ry = 3.1 (WD model) and that from 
Type-II supernovae (SN model). The arrow at ~ 9.43 fim indicates a 2<j 
upper limit. Jl 148+5251 was observed with MIPS at 70 fim but with no 
detection. The arrow indicates a 2cr upper limit. The current observations 
cannot distinguish between these two dust models. 



agree well at wavelengths A > 10 fim, but differ by a factor of 
a couple in the optical /NIR bands 0.1- 10 /im. It appears that 
WD model would require a higher dust-to-gas ratio than that 
of the Milky Way by a factor of a few in order to produce the 
observed SED of Jl 148+5251. In such case, the WD model 
would produce a stronger peak at the 9.7 fim silicate feature 
by a similar factor relative to the SN model. Observation of 
J 1148+5251 at ~ 9.1 fim would be helpful in distinguishing 
between these two models. However, the current data-point 
at that wavelength (observed by MIPS at 70 fim) is only a 2a 
upper limit, insufficie nt to c o nstrain the model. 

As discussed in § 12.2.41 Elvis et all (120021) suggest that 
quasars may also be copious producers of dust, as condensa- 
tion in quasar outflows is similar to dust formation in the en- 
velopes of AGB stars. In this case, the dust extinction curve 
would be similar to the WD model. However, the contribu- 
tion of dust from quasar winds i n the early Universe remains 
unknown. The recent report by IStratta et al.1 (120071) of dust 
extinction in the host galaxy of GRB 050904 at z = 6.3 in- 
dependently supports a supernova-origin dust model. If both 
supernovae from starbursts and quasar outflows play impor- 
tant roles in dust production in these young quasar systems at 
z ~ 6, then the resulting dust extinction curve would be inter- 
mediate to the WD and SN models shown in Figure [TT] More 
observations and deeper surveys for dust in high-z galaxies 
and quasars will be necessary to constrain dust formation 
mechanisms at early cosmic times. 

4.5. Dust-to-gas Ratios 

It has been suggested that ULIRGs have dust-to-gas ratio 
(mostly col d gas as in our m odelin g) close to th at of the Milky 
Way (e.g.. iDunne & Ealesl l200lt iKlaas et all 120011) . While 
some observations of obscurred, X-ray selected AGNs seem 
to suggest that the dust-to-gas ratio in these obj ects has a wide 
range , from ~ 10~ 3 to afewof MW value (e.g..lMaiolino et al.l 
I200U iKuraszkiewicz et al.ll200l lHall et al.ll2006l) . Although 
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FIG. 12. — A parameter study of the dust-to-gas ratio relative to Milky 
Way value, for the Cold Phase Gas (CPG, top panel) and Hot Phase Gas 
(HPG, bottom panel), respectively. This figure shows that the emergent SED 
is more sensitive to the dust-to-gas ratio of the HPG than to the CPG owing 
to the large covering factor of the HPG. In the regular RT calculations, the 
default values are 1 and 1 % Milky Way value for CPG and HPG, respectively. 



this range may apply mainly to the hot, fully or partially ion- 
ized gas in the circumnuclear regions of the AGNs, it would 
be interesting to see how such a range affects the output SED 
in our calculations. 

Figure[T2]shows the comparsion of the emergent SEDs with 
different dust-to-gas ratios in a wide range, for both cold (top 
panel) and hot-phase gases (bottom panel), respectively. The 
values in the plot are relative to the MW value. A change 
of two orders of magnitude in the dust-to-gas ratio of the cold 
phase gas results in a difference in the cold dust bump only by 
a factor of a few. However, a similar change in the dust-to-gas 
ratio in the hot phase gas would result in substantial differ- 
ence in the output SED. This study shows the extinction is 
dominated by the hot dust, which has a much larger covering 
factor that the cold dust. We find that using 1 and 1 % of MW 
value for the cold and hot phase gas, respectively, reproduce 
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the observation of SDSS Jl 148 reasonably well. Therefore, 
we elect to use them as default values for our standard RT 
calculations. Note this choice of 1% of MW value for the hot, 
diffuse gas has little effect for the starburst or ULIRG phase, 
because a substantial fraction of gas during that stage is in the 
cold clouds. However, it may affect the SEDs of major quasar 
phase, as the gas is heated by the AGN. We should point out 
that in our IS M model, we do no t have the warm phase as in 
the picture of iMcKee & Ostrikerl dl977l) . Such a warm phase 
medium may have a dust-to-gas ratio similar to that of MW. 
We will explore such a treatment in future work. 

5. THE DUST DISTRIBUTION 

The distribution of dust is essential to investigating dust for- 
mation mechanisms and heating sources. Figure Q~3] shows 
the projected spatial distribution of the dust density and tem- 
perature in the system at z = 6.5. The projected quantity is 
calculated with J f(l)dl/L, where /(/) is density or temper- 
ature distribution along the line of sight, while L is the total 
length of the sight line. Gravitational torques in the interac- 
tion produce strong shocks and tidal features, making the dust 
distribution highly inhomogeneous and clumpy. The dynamic 
range in density can be up to six orders of magnitude. Both 
the dust density and temperature peak in the central regions 
near the AGN. 

The detailed radial distributions of the dust mass and tem- 
perature are quantified in Figure [14] The blue and red curves 
in this figure represent CPG- and HPG-dust, dust associ- 
ated with cold- and hot-phase gas, respectively, as defined 
in § 12.2.31 The CPG-dust has high density but a small vol- 
ume filling factor; it is generally optically thick to radiation. 
On the other hand, the HPG-dust has a lower density but fills 
> 99% of the volume, and is generally optically thin. Both 
the dust mass and temperature vary with distance from the 
central AGN. The cold dust is typically surrounded by hot 
dust. In the inner hundred parsecs, the cold cores are highly 
condensed and are optically thick even to the hard radiation 
from the black hole. However, the HPG-dust at the edges 
surrounding these cold cores is heated directly by the central 
AGN. This is indicated by the power-law temperature profile 
in the lower panel of Figure [14] T oc R~ l l 2 as expected from 
Equation ©: T 4 (R) oc E(R) oc R~ 2 . The hottest dust is heated 
up to ~ 1200 K, which is below the dust sublimation temper- 
ature - 1600 K dHonig et al.ll2006l) . Cooler dust (T ~ tens to 
hundreds K) is distributed in an extended region from ~ 100 
pc to several kiloparsecs. The heating sources come both from 
stars and the central AGN. Beyond 10 kpc, the dust tempera- 
ture drops to below 100 K. Note that the minimum tempera- 
ture in Figure [14] is ~ 20 K, which is the CMB temperature at 
z=6.5. 

Figure Q3] shows histograms of the dust mass and temper- 
ature. The dust temperature ranges from ~ 10- 10 3 K. The 
amount of the hottest dust (~ 10 3 K) is — 1 .4 x 10 2 M Q , while 
that of the cold dust (T < 10 2 K) is - 1.4 x 10 8 M Q . This 
is in agreement with the estimates of the dust detected in the 
host of J 1148+52 5 1 dBertoldi et alj2003at iBeelen et al.ll2006t 
Uiang et al.ll2006l) . 

The large amount of cold dust (~ 1.4 x 10 8 M Q ) lo- 
cated within 3 kpc from the AGN provide s efficient cool- 
ing fo r the formation of molecular gas. In Narayana n et al.l 
(2007), we calculate carbon monoxide emission using a 
non-local thermodynam ic equilibrium radiative transfer code 
(Naravanan et al. 2006a b), and find that CO gas forms in this 
region, with a total mass of ~ 10 10 M Q , similar to observa- 




-40 -20 20 40 60 
X (kpc) 



FIG. 13. — Maps of the projected density (top panel) and temperature 
(bottom panel) of the HPG-dust (dust associated with hot phase gas). The 
origin of the map is the location of the central quasar, and the coordinates 
are comoving. Note that the gas only determines the distribution and mass of 
the dust; the dust temperature is not associated with the gas t emper ature but 
is calculated self-consistently from the radiation field (see j) l2.2.3l for more 
details). 



tions bv l Walter et al.l d2004l) . This suggests that cold dust and 
molecular gas are closely associated, and both depend on the 
star formation history. Our results show that significant metal 
enrichment takes place early in the quasar host, as a result of 
strong star formation in the progenitors, and that intense star- 
bursts (SFR > 10 3 M Q ,yr _1 ) within < 10 8 yr are necessary to 
produce the observed properties of dust and molecular gas in 
the host of J 1 148+5251, a conclusion which is also supported 
by the analytical models of lDwek et alj d2007l) . 
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FIG. 14. — Radial distribution of dust mass (top panel) and temperature 
(bottom panel) as a function of the distance from the central AGN. The blue 
and red curves represent the CPG- and HPG-dust, respectively, while the 
black curve in the bottom panel is a power-law temperature profile, T oc R~ l l 2 
as expected from heating by the central AGN. 
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FIG. 15. — Mass-temperature distribution of both the CPG- (in blue) and 
HPG-dust (in red). In the bottom panel, M c m, M warm , and Mh t represent the 
total mass of cold-, warm- and hot-dust with T < 100 K, 100 < T < 1000 K, 
and 1000 < T < 1200 K, respectively. 
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FIG. 16. — Evolution of the SEDs of the quasar system and its galaxy 
progenitors in the observed frame. The colored curves represent S EDs from 
z ~ 14 to z ~ 5.2, while the black dots are again observations from J iang et alj 
12006), as described in the legend. Absorption of the Lyman line series and 
continuum by the intergalactic medium iMadau 1995) is taken into account, 
which results in a sharp drop at short wavelengths. 



6. EVOLUTION OF THE QUASAR SYSTEM 



6.1 



. Transition from cold to warm ULIRG 

(l2007h . we show that the host of the z ~ 6 quasar 

14, and the 



In lLi et alJ 

undergoes hierarchical mergers starting from z ■ 
quasar descends from starburst galaxies. The evolution of the 
SEDs of the system in the observed frame is shown in Fig- 
ure [16] Note that absorption of the Lyman lines and con- 
tinuum by the intergalactic medium (Madaui ri995h is taken 
into account. During early stages at z ~ 14, the SED of the 
quasar progenitor shows only a cold dust bump that peaks 
aro und ~ 60/xni, which is ch aracteristic of starburst galax- 
ies (Sanders & Mirabel 1996). As the system evolves from 
starburst to quasar phases, dust extinction in the UV-optical 
bands increases, boosting the emission reprocessed by dust at 
wavelengths longward of 1 /im. The gradually increasing ra- 
diation from the accreting black holes heats the nearby dust 
to high temperatures, contributing to the hot dust bump which 
peaks around ~ 3/j.m (rest frame). At the maximum quasar 
phase at z — 6.5 indicated by the red curve, the hot dust bump 
SED reaches its peak with a temperature of ~ 1200 K, as we 
have seen in the previous section. Such an SED represents 
lu minous, blue quasars in the samples of Jiang et alj d2006l) 
and lRichards et al.l J2006). As the system ages and reddens in 
the post-quasar phase (indicated by the black curve), the total 
luminosity of the system drops. However, there is still some 
residual hot dust, and the infrared luminosity is dominated 
by the NIR and MIR. This resembles the class of infrared- 
bright, optically-r ed quasars found in recent surveys (e.g., 
iBrand et al.ll2006l) . Overall, the evolution of the SED from 
a starburst to a quasar can be characterized by the slope of the 
infrared SED (3 - 50 /im), as it decreases from the starburst to 
quasar phases owing to the increase of NIR emission from the 
hot dust heated by the AGN. 
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FIG. 17. — Evolution of rest-frame infrared luminosities of the quasar sys- 
tem (top panel) and the rest-frame color fis^m/ fmpt,m (bottom panel). In the 
top panel, the solid curves represent the luminosity in near-IR (1 — 10/im), 
mid-IR (10-40 /im), far-IR (40-120 /im), and IR (8-1000 /im), respec- 
tively. The grey shades indicate regimes for LIRG (Lir > 10 u Lq), ULIRG 
(Li r > 10 12 L W ), and HLIR G (Lir > 1O 13 L ), respectively, as classified 
by [Sanders & Mirabel 1 1996). The ratio /25(jm//60/jm in the bottom panel 
is an indicator of the co ldness of the SED with a critical value of 0.3 
(Sanders & Mirabel 1996). The quasar system in our model evolves from 
"cold" ULIRG (/25 M m//60 M m < 0.3) to "warm" ULIRG (including HLIRG, 
flSfim/ f 60 iim > 0.3) as it transforms from starburst to quasar phases. 



Infrared luminosities are powerful tools to study starburst 
galaxies and quasars. From the results in the previous sec- 
tions, we see that the emission in near-IR (1 - 10 /im), mid-IR 
(10-40 /im) and far-IR (40- 120 /im) comes from re-emission 
by hot, warm, and cold dust, respectively. Note that the mean- 
ing of FIR varie s in the li t eratur e. Here we use the definition 
of FIR given by ICondonl (Il992|) in order to comp are our re- 
sults with the observations by ICarilli et"aTI (120041) who used 
the same FIR range. 

Using L IR (8 - 1000 ^m). ISanders & Mirabell dl996b classi- 
fied infrared luminous galaxies into three categories, namely 
luminous infrared galaxy (LIRG, Lir > 10 L ), ultra-LIRG 
(ULIRG, Lir > 1O 12 L ), and hyper-LIRG (HLIRG, Lir > 



10 13 L Q ). Figure [T7] (top) shows the evolution of the infrared 
luminosities of the quasar system in our simulations. The 
luminosities increase with the star formation rate and black 
hole accretion rate. The system is ultraluminous most of the 
time, with periods in HLIRG phases associated with bursts 
of star formation or quasar activity. When the last major 
mergers take place between z ~ 9-7.5, the strong shocks and 
highly concentrated gas fuel rapid star formation and black 
hole growth. This also produces a large amount of dust heated 
by the central AGN and stars. As a result, the infrared lumi- 
nosities increase dramatically, pushing the system to HLIRG 
class. 

The infrared luminosities strongly peak at z ~ 8.7 when the 
star formation rate reaches > 10 4 M Q yr _1 , suggesting a sig- 
nificant contribution from stars in heating the dust to emit 
in the range 1 - 1000/im. During the major quasar phase 
(z ~ 7.5-6) the star formation declines to ~ 10 2 M Q yr" 1 , and 
the emission in NIR and MIR outshines that in FIR, demon- 
strating that AGN can play a dominant role in heating the 
dust and producing NIR and MIR emission. Furthermore, the 
AGN can also contribute to FIR emission. For example, while 
the star formation rate drops by a factor of ~ 500 from z ~ 8.7 
(when star formation rate peaks) to z ~ 6.5 (when black hole 
accretion rate peaks, see FigureQ]), the Lfir declines by only a 
factor of ~ 50, indicating substantial contribution to the Lfir 
by the AGN. This significant AGN contribution has important 
implications for estimating the star formation rate. We will 
discuss this issue at length in the next section. During the 
"post-quasar" phase at z < 6, as a result of feedback which 
suppresses both star formation and black hole accretion, and 
gas depletion which reduces the amount of dust, the infrared 
luminosities drop rapidly. 

The flux ratio /2s Mm / '/isO/un is a color indicator for the cold- 
ness of the SED, and c an be used as a diagnostic for AGN ac- 
tivity, as suggested by Sanders & Mirabel ( 1996). For exam- 
ple, starburst galaxies usually have cold colors /25 Mm //60^m < 
0.3, while quasars are warm with fzSfim/ feOfim ^ 0.3. The 
evolution of the color fiSuml fobum of the simulated quasar 
system is shown in Figure [T7] (bottom panel). The color in- 
dex is below 0.3 most of the time, i.e. during both starburst- 
and post-quasar phases. However, it rises above 0.3 dur- 
ing the short quasar phase. This figure clearly demonstrates 
that the system evolves from a cold to warm ULIRG (in- 
cluding HLIRG) as it transforms from starburs t to quasar 
phases. Similar trend is also reported by Chakrabar ti et alj 
d2007al) . Our results provide further theoretical support for 
the st arburst-to-quasar conjecture, as su ggested by observa- 
tions (iindiriAMlr*il[l99l HcovHSiS)!). 

6.2. AGN Contamination and the SFR - Lfir Relation 

AGN contamination in infrared observations of dusty, star- 
forming quasar systems has been a long-standing problem. In 
our simulations, we see that the contributions from AGNs and 
stars both vary according to the activity of these two popula- 
tions. As shown in Figure [18] during the starburst phase, not 
surprisingly nearly all the infrared emission comes from dust 
heated by stars. However, during the peak quasar phase, the 
contribution from the AGN dominates the infrared light pro- 
duction in the system. This result has significant implications 
for the interpretation of observational data, such as estimates 
of the star formation rate. 

In particular, owing to the lack of other indicators such 
as UV flux or Ha emission, in observations of high-redshift 
objects the far-infrared luminosity is commonly used to es- 
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FIG. 18. — Contribution to the infrared luminosities, Lnir, Lmir, and Lfir, 
from both AGN and stars, respectively. During the starburst phase, stars are 
the main heating source for the infrared emission. However, during the peak 
quasar phase, AGN heating dominates. 



timate the star formation rate, assuming most or all of Lfir 
is contributed by young stars. For example, Jl 148+5251 has 
L FIR ~2x 1O 13 L , and t he star formation rate is estimated 
to be ~ 3 x K^Moyr" 1 dBertoldi et alj|2003at ICarilli et al.1 
2004). However, these studies cannot rule out a substantial 
contribution (e.g. ~ 50%) to Lfir by the AGN. If Lfir is in- 
deed heavily contaminated by AGN activity, then the SFR - 
Lfir relation would be invalid, and the actual star formation 
rate would be much lower. 

To demonstrate this, we show in Figure [19] the SFR - Lfir 
relation in three cases: (1) only black holes are included as 
a heating source of the dust (top panel); (2) only stars are 
included (middle panel); and (3) both black holes and stars are 
included (bottom panel). The top panel shows that there is no 
correlation between the SFR and the Lfir which is produced 
by AGNs. However, if all the Lfir is produced solely from 
stars, then there is a tight, linear correlation. If the Lfir is 
contributed by both AGN and stars, then the correlation would 
change. The bottom panel in Figure [19] shows that during the 
starburst phase where stars dominate dust heating, there is still 
a correlation between SFR and Lfir with modified slope and 
normalization. However, during the peak quasar phase where 
the AGN dominates dust heating, there is no correlation at 
all. To summarize, Figure [l9]gives the following SFR - Lfir 
relations: 



SFR 



M Q yr 
SFR 



(L \ 101 

T = 2.3x 10- 10 f j^j (stars-only), (21) 



M Q yr- 



0.76 



r =4.9 x 10"' (AGN+ stars). (22) 



Here, the FIR is in the range of 40- 120/i m (rest frame). 
This correlation in the stars-only case is very close to that 
foundby lKelmicurl(ll998blK ^^-r ~ 1.7 x 10 -10 ^ for star- 
burst galaxies. Note the slight difference in the normaliza- 
tion might owe to di fferent stellar I MFs us ed - we us e the top 
heavy Kroupa IMF dKroupall2002|), w hile iKennicutt] (Il998bl) 
adopts a Salpeter IMF (Salpetei IT955h . Also, the Lpm in that 
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FIG. 19. — Evolution of the SFR - Lfir relation in our simulations. Here 
we consider the relation in three cases where Lfir is contributed by AGNs 
only (case 1, top panel), stars only (case 2, middle panel), and both AGNs and 
stars (case 3, bottom panel). The colored filled symbols indicate the system at 
different redshifts, while the black open diamond represents the model quasar 
at z = 6.5. The black curve in the middle panel is the least-squares fit to all the 
data, while that in the bottom panel is the fit to data points from z ~ 14—7.5 
only (starburst phase). Case 1 has no SFR - Lfir correlation; case 2 has a 
tight, linear correlation similar to that used in observations; and case 3 has a 
non-linear correlation. 
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FIG. 20. — Evolution of the dust mass and temperature as a function of 
redshift. The dust is heated by both stars and AGN. As the system evolves 
from starburst to quasar phases, the amount of hot dust increases. The dust 
reaches the highest temperature during the peak quasar phase. After that, the 
dust cools owing to the decline of AGN radiation. 



work refers the infrared luminosity integrated over the range 
of 8-1000/xm. 

In our simulation, the model quasar at z = 6.5 has Lfir ~ 
1.8 x 1O 13 L , which is c lose to that of Jl 148+5251 derived 
from radio observation dBertoldi et alj l2003at ICarilli etail 
2004). If we assume a perfect SFR — Lfir correlation as 
in the middle panel of Figure [19] (Equation lF2TI ). then the 
S FR would be ~ 4 x 1 3 M f 7)yr~ 1 , similar to the estimate 
of IBertoldi etail (l2003ah : ICarilli et al.1 (l2004h . However, the 
SFR of the model quasar from our simulations is about ~ 
112MQyr _1 , which is significantly lower than that derived 
from an ideal SFR - Lfir correlation. This discrepancy is 
largely due to significant contribution to the Lfi r by AGN ac- 
tivity Sev eral studies bv IBertoldi et alJ (l2003ah : ICarilli etalJ 
(2004) and lWang et al.l d2007l) find that most z ~ 6 quasars do 
not seem to show a Lfir - Lb correlation, but instead fol- 
low a strong Lfir - Lradin correlatio n as measured in local 
star-forming galaxies (Condonl ll992l) . These authors there- 
fore suggest that these z ~ 6 quasars are strong starbursts, and 
that most of the Lfir may come from stars. However, the 
Lfir - L ra di correlation may not guarantee the stellar origin 
of the Lfir as the physical basis for this correlation is un- 
known. Our results suggest that in quasar systems, Lfir alone 
may no longer be a reliable estimator for star formation rates; 
other diagnostics should also be considered. We will study 
this topic with more detail in a future paper by investigating 
the correlations in a multiple luminosity plane, Lb - Lfir - 
Lx-ray - L la dio, of starburst galaxies and quasars. 

6.3. Evolution of the Dust Properties 

Figure |20] shows the evolution of the dust properties, in- 
cluding the dust mass and temperature, as a function of red- 
shift. As the system evolves from starburst to quasar phases, 
the amount of hot dust increases accordingly owing to the en- 
hanced heating from the central AGN, as well as replenish- 
ment of gas /dust from incoming new galaxy progenitors dur- 
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FIG. 21. — Relations between rest-frame flux fi^ m and hot dust mass 
A^hot (t°P panel), and /sopm) an d c °ld dust mass M co i,j (bottom panel). The 
black curve is the least-squares fit to all the data in top panel, and data during 
starburst phase at redshift 7, ~ 14—7.5 in bottom panel. 



ing mergers. The amount of cold dust (T < 100 K) reaches a 
maximum at z ~ 8.8 when the last major merger takes place, 
then decreases steadily as the dust is heated to higher tem- 
peratures, or as the gas is consumed by star formation and 
black hole accretion. Strong star formation is able to heat 
the dust nearly to T ~ 1000 K. However, the hottest dust at 
T ~ 1200 K is associated only with the peak quasar phase at 
z ~ 6.5 when the hot dust is directly heated by the AGN. After 
that, the amount of hot dust drops dramatically as a result of 
the rapidly declining radiation from the central AGN. 

In observations, the dust mass is usually determined by 
assuming that the dust radiates as a black- or grey -body 
dHughes et al.lll997l;lJiang et al.ll2006l) . In the emergent SEDs 
from our radiative transfer calculations, the cold dust bumps 
at different redshifts appear to peak around ~ 50/im (rest- 
frame), while the hot dust bumps peak around ~ 3/im. Fig- 
ure [2U shows the relations between the dust mass and fluxes 
at these two wavelengths. The top panel shows that M^ot 
(T > 10 3 K) increases with/ 3Aim flux, while M coU (T < 10 2 K) 
correlates with flux linearly during the starburst phase: 
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FIG. 22. — Evolution of the relation between rest-frame 3.5 (im luminosity 
and B-Band luminosity (4400 A) for the quasar system in our simulations. 
The filled circles represent starburst phases, squares represents quasar phases, 
and triangles represent post-quasar phases. Colors indicates redshift of the 
object. 
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However, during the quasar phase, there appears to be an 
excess in fso^m, so the M co id - correlation does not ap- 
ply. These results suggest that rest-frame 3/im and 50/im 
might serve as good diagnostics for dust, the former for hot 
dust in quasar systems, and the latter for cold dust in starburst 
galaxies. 

It has been suggested bv lJiang et al.l (|2006) that the NIR-to- 
optical flux ratio may be used to probe dust properties. From 
a sample of thirteen z ~ 6 quasars observed with Spitzer, these 
authors find that two quasars have a remarkably low flux ra- 
tio (rest-frame 3.5^m to B-band) compared to other quasars 
at different redshifts, and they are also weak in FIR. Further- 
more, such a low flux ratio was not seen in low-redshift quasar 
counterparts. These findings lead to the suggestion that these 
two quasars may have different dust properties from others 
(iJiang et al.ll2006h . 

From our model, we find that the dust properties are associ- 
ated with the evolutionary stages of the host. Figurel22lshows 
the evolution of the ratio Lnir/Lb from the quasar system at 
different stages of its life. Early on, during the starburst phase, 
the ratio Lnir/Lb is low. As the system proceeds to its peak 
quasar activity, the heating from the central AGN increases 
the hot dust emission around 3 /im, boosting the Lnir/Lb ra- 
tio. After that, in the post-quasar phase, as the radiation from 
the AGN declines, the hot dust emission drops rapidly, re- 
sulting in a lower Lnir/Lb- According to our model, there 
might be two po ssible explanations for the two outliers in the 
Jiang et al. (2006) sample. One possibility is that they may 
be young quasars that are still in the starburst phase but have 
not yet reached peak quasar activity, so the light from star for- 
mation may be dominant or comparable to that from the ac- 
creting SMBH still buried in dense gas. This may explain the 



low NIR flux, as well as the B-band luminosity and the nar- 
row Lya emission line, which are primarily produced by the 
starburst. From Figure[l7] we note that there are "valleys" be- 
tween major starbursts (or mergers) where the infrared lumi- 
nosities are low. This may explain the weak FIR in these two 
quasars. Another possibility is that these two objects may be 
in the post-quasar phase, which would also have low Lnir/Lb 
and FIR (see also Fi gure [T71 and \T% . 

However, our model predicts that the same evolution should 
apply to every luminous quasar formed in gas-rich mergers. It 
cannot explain statistically why such a low Lnir/Lb ratio is 
only found in the z ~ 6 quasars, but not at low redshifts. Per- 
haps these two quasars indeed have unusual dust properties 
that cannot be explained by our current model. More observa- 
tions of dust in high-z objects are n eeded to resolve th e issues 
surrounding the IR-weak quasars in lJiang et al.l (120061) . and to 
test our model. 

7. DISCUSSION 

Our three-dimensional Monte Carlo radiative transfer code 
ART 2 makes it possible to calculate self-consistently the ra- 
diative transfer and dust emission from galaxies and quasars 
which have large dynamic ranges and irregular geometries. 
It works efficiently on hydrodynamics simulations of galax- 
ies and mergers, and has flexible grid resolution that can be 
adjusted to any SPH resolution in hydrodynamic simulations. 

To date, only a limited number of RT calculations 
have been performed o n galaxy mergers dJonssonl 2006; 
IChakrabarti et al.ll2007allbl). In particula r, the parallel SUN- 
RISE code developed by iJonssonl (2006) marked a milestone 
in this direction. It has an adaptive grid similar to ours and 
works efficiently on the arbitrary geometry in galaxy merg- 
ers. However, radiative equilibrium is not yet included in this 
code, so SUNRISE can not currently calculate dust emission 
self-consistently. 

On the other hand, IChakrabarti et all d2007allbh had a sim- 
ilar app roach to ours. They apply the code of Iwhitney et alJ 
(l2003bl) . w hich employs the same ra diative equilibrium algo- 
rithm as in Biorkman & Wood d2001l) . to hydrodynamic sim- 
ulations of galaxy mergers with black holes. However, the 
spherical, logarithmically-spaced grid used by these authors 
is not optimal for describing an inhomogeneous distribution 
of gas and dust with multiple density centers in galaxies and 
mergers, as exemplified by Figure [3] Moreover, our method- 
ology is based on a treatment of the multiphase ISM that in- 
cludes extinction from dust in the diffuse phase, which was 
ignored in the work of IChakrabarti et all (|2007a b). Further- 
more, we adopt an AGN input spectrum that includes emis- 
sion fr om dust on scales that the simulations do not resolve, 
unlike Chakrab arti et al.l d2007albl) who employ a power-law 
AGN spectrum. The tests described above demonstrate that 
the near- and mid-IR emission is sensitive to the resolution of 
the RT calculation (Figure|7]i, extinction by dust in the diffuse 
phase of the ISM (Figure[8]), and the templat e spectrum of the 
AGN (F igurelTOll. The SEDs in the work of lChakrabarti et alJ 
d2007allbl) show only cold dust emission characteristic of a 
starburst during all evolut ionary phases, even whe n the quasar 
peaks (e.g., Figure 20 in Chakrabart i et al.l2007bl) . and, in par- 
ticular, do not at any time exhibit hot dust emission. We at- 
tri bute the difference s between our computed SEDs and those 
of Chakra barti et alJ d2007albh mainly to the different grid 
method employed, the handling of extinction by the diffuse 
ISM, and our choice of an AGN input spectrum that includes 
emission on scales that cannot be resolved in our hydrody- 
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namic simulations. 

In our radiative transfer calculations, we do not explicitly 
include a dust torus near the AGN. The resolution of our 
RT grid is constrained by that of the hydrodynamic simula- 
tions, which is of ~ 30 pc. As a result, we can not resolve 
the ne ar vicinity of the AGN. In the Unification Scheme for 
AGN dMiller & Antonucci| [l983). a parsec-scale dust torus, 
both optically and geometrically thick, has been used to suc- 
cessfully explain the different appearances of Type 1 and Type 
2 AGNs. For example, Type 1 AGNs are typically viewed 
face-on and show blue UV bump in their spectra, while Type 
2 AGNs are viewed edge-on without the blue UV bump. 
It is believed that the NIR and MIR emission of an AGN 
comes from the hot dust in the torus structure, as suggested 
by various observations, for examples, the MIR and NIR de- 
tection of a Seyfert 2 galaxy NGC 1068 (iJaffe et alJ 12001 
IWittkowski et alj | 2004l), and the MIR emi ssion of a lensed 
quasar Q2237+0305 by Wvith eet al.1 d2002l) in the CfA Red- 
shift Survey (Huchra et al. 1985b- The radiative transfer mod- 
eling by Honig et al.] (l2006h . which concentrates on the in- 
ner parsecs of the AGN, has successfully reproduced the SED 
of NGC 1068 in the NIR and MIR with a three-dimensional, 
clumpy dust torus. Such a scale is, however, below the resolu- 
tion of our modeling, which follows the dynamical evolution 
of the quasar system on a galactic scale. Moreover, it is not 
clear whether such a well-defined torus structure would exist 
in luminous quasars at high redshift s. Nevertheless, we ha ve 
adopted a composite AGN spectrum (Hopkin s et alj2007dl) as 
a sub-resolution recipe for the IR emission from the hottest 
dust near the AGN. Such an approach is supported by the 
fact that the composite spectrum represents th e average near- 
to mi d-IR emission in thousands of quasars (Richar ds et alJ 
2006), and that the torus does not contribute much to the far- 
IR emission, which comes mainly from cold dust on kilopar- 
sec scales. Therefore, for spatial scales of interest, the results 
should be reliable, and they are expected be similar to those 
from calculations with a torus included. We defer such a cal- 
culation to the future when sub-parsec resolution is feasible 
in hydrodynamic simulations of quasar formation. 

We should point out that dust destruction is not ex plicitly 
included in our modeling. It ha s been suggested (e.g.,[M c Kee 
1989; Draine&McKeei l 19931: iDraind 120031; iNozawa et all 
2006) that dust destruction might be efficient in non-radiative 



shocks by non-thermal and thermal sputtering owing to a 
high shock velocity (> 100 km s" 1 ) and a high gas temper- 
ature (> 10 5 K). However, the sputtering timescale for this 
temperature is of ~ 3 x 10 6 yr for a gas density n = 1 cm" 1 , 
and even shorter if the gas density is higher dBurke & Silkl 
1 19741; [Reynolds et al.l 1 1997b . It has also been suggested 
that dust grains may b e destroyed behind radiative shocks 
dTodini &Ferraral 120011: iBianchi & Schneider! 120071) . How- 
ever, it is expected that the efficiency cannot be significant at 
high r edshifts as magnetic fields are likely weak ( Gne din et alJ 
2000). In our dust treatment in the RT analysis, the dust 
is re-calculated based on the gas content in every snapshot 
from the hydrodynamic simulations, which has a time inter- 
val of 2 Myr. This is equivalent to dust destruction within this 
timescale. Moreover, we assume different dust-to-gas ratios 
for the cold and hot gas, which accounts for dust survival in 
these different phases of the interstellar medium. 

One limitation of our ART 2 code is that it does not in- 
clude the transient heating of small dust grains (size < 
200A), which may cause temperature fluctuation and pro- 



duce enhanced near-infrared emission dLi & Draind l200lt 
iMisseltet ail 120011) . in particular for the PAH feature. Al- 
though PAH is not studied in this present work, the code 
can be improved in the future by implementing a temperature 
distribution for those small grains. Moreover, one potential 
caveat of our RT work is that it is a post-processing procedure. 
It would be ideal to couple the radiation with hydrodynamics; 
i.e., having the radiative transfer done simultaneously with the 
hydrodynamic simulations (as in e.g., lYoshida et al.ll2007h as 
the dynamical evolution of the gas might affect the formation 
and properties of the dust. However, this effect is currently 
uncertain because such a technique is not yet available for 
dust emission in galaxy simulations. Although ART 2 is not 
parallelized, it works fairly efficiently owing to the adaptive 
refinement tree. 

We have done extensive resolution studies which include 
both the number of photon packets and the refinement level 
of the grid, and parameters for the ISM model such as the 
cloud mass spectrum and size distribution, as well as differ- 
ent dust models. We find convergence of the modeled SEDs 
with parameters within observational ranges (see § 2 and § 3). 
Therefore, we conclude that the RT calculations we present 
here with ART 2 are robust, and that our model is more self- 
consistent, more realistic, and more versatile than previous 
approaches. 

8. SUMMARY 

We have implemented a three-dimensional Monte Carlo 
radiative transfer code, ART 2 - All-wavelength Radiative 
Transfer with Adaptive Refinement Tree, and use it to cal- 
culate the dust emission and multi-wavelength properties of 
quasars and galaxies. ART 2 includes the following essential 
implementations : 

1 . A radiative equil i brium algorithm developed by 
iBiorkman & Wood! d2001l) . It conserves the total pho- 
ton energy, and corrects the dust temperature without 
iteration. This algorithm calculates dust emission effi- 
ciently and self-consistently. 

2. An adaptive grid s cheme in 3-D C artesian coordinates 
similar to that of Jonsson (2006), which handles an 
arbitrary geometry and covers a large dynamic range 
over several orders of magnitude. This is indispensable 
for capturing the inhomogeneous and clumpy density 
distribution in galaxies and galaxy mergers. It easily 
achieves a ~10 pc-scale grid resolution equivalent to a 
~ 10000 3 uniform grid, which is prohibitive with cur- 
rent computation schemes. 

3. A two-phase ISM model in which the cold, dense 
clouds are embed ded in a hot, diffuse medium in pres- 
sure equilibrium (Springel & Hernquist 2003a). More- 
over, the cold clouds follo w a mass spec trum and a size 
distribution similar to the iLarsonl d 198 lh scaling rela- 
tions inferred for giant molecular clouds. This model 
ensures an appropriate sub-grid recipe for the ISM 
physics, which is important for studying dust proper- 
ties in galaxies. 

4. A supernova-origin dust model in which the dust is 
produced by Type-II supernovae, and the size distribu- 
tion of grains follows that derived by Todini & F erraral 
d2001l) . This model may be especially relevant for dust 
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in high-redshift, young objects. In traditional dust mod- 
els, the dust is produced by old, low mass stars such as 
AGB stars, which are typically over 1 Gyr old. How- 
ever, quasar systems at ~ 6.4 are only a few hundred 
Myr old, and there would be insufficient AGB stars to 
produce the abundant dust as observed in J 1 148+5251. 

5. The input spectra include those from both stars, cal- 
culated using STARBURST 99 dLeitherer et al.1 119991: 
Vazquez & Leitherer] 120051) . and black holes repre- 
sented by a composit e AGN template developed by 
Hopk ins et al.l (l2007dl) . The composite spectrum in- 
cludes a broken power-law as the intrinsic black hole 
spectrum, and infrared components that come from the 
torus near the AGN. The torus structure is unresolved 
in our hydrodynamic simulations. This template is a 
sub-grid recipe to include the hot dust emission within 
pc-scales which is below the resolution of our hydrody- 
namic simulations. 

ART 2 works efficiently on hydrodynamic simul ations of 
galax ies and quasars performed with GADGET2 (Springel 
2005). By ap plying ART 2 to the quasar simulations of 
Li et alJ ([2007), in which luminous quasars at z > 6 form 
rapidly through hierarchical mergers of gas-rich galaxies, we 
are able to calculate the multi-wavelength SEDs (from opti- 
cal to millimeter) and dust properties of the model quasar at 
z ~ 6.5 and its galaxy progenitors at even higher redshifts. 

We find that a supernova-origin dust model may be able to 
explain the dust properties as observed in the high-redshift 
quasars. Our calculations reproduce the observed SED 
and properties such as the dust mass and temperature of 
J 1 148+5251, the most distant Sloan quasar. The dust and in- 
frared emission in quasar hosts are closely associated with 
the formation and evolution of the system. As the system 
transforms from starburst to quasar phases, the evolution of 
the SEDs is characterized by a transition from a cold to warm 
ULIRG. During the starburst phase at z > 7.5, the SEDs of the 
quasar progenitors exhibit cold dust bumps (peak at ~ 50 /im 
rest frame) that are characteristic of starburst galaxies. Dur- 
ing the peak quasar phase (z ~ 7.5-6), the SEDs show a 
prominent hot dust bump (peaks at ~ 3 /im rest frame) as ob- 
served in luminous quasars at both low and high redshifts (e.g, 
iRichards et afl200ftl.Tianp et al.ll2006h . 

Furthermore, we find that during the quasar phase, AGN ac- 
tivity dominates the heating of dust and the resulting infrared 
luminosities. This has significant implications for the inter- 

2 http://www-star.st-and.ac.uk/~kw25/research/montecarlo/montecarlo.html 

3 http://gemelli.colorado.edu/~bwhitney/codes/codes.html 

4 http://www.ucolick.org/~patrik/sunrise/ 



pretation of observables from the quasar host. The hottest dust 
(r > 10 3 K) is present only during the peak quasar activity, 
and correlates strongly with the near-IR flux. The SFR - Lfir 
correlation depends sensitively on the relative heating con- 
tributed from AGN and stars. If the dust heating is dominated 
by stars, as in the starburst phase, then there is a tight, linear 
SFR - Lfir correlati on similar to the one widely used to inter- 
pret observations dKennicuttlll998bl) . However, if both AGN 
and stars contribute to the dust emission in the far-IR, then 
the correlation becomes non-linear with a modified normal- 
ization. If the AGN dominates the Lfir, as in the peak quasar 
phase, we find no correlation at all. Finally, we find correla- 
tions between dust masses and rest-frame fluxes at 3 //m and 
50 /^m. The fc^m flux may serve as a good diagnostic for 
hot dust in quasars, while /so^m may be used to estimate the 
amount of cold dust in starburst galaxies, as /so^m correlates 
linearly with the cold dust mass. 

Our model demonstrates that massive star formation at 
higher redshifts (z > 10) is necessary in order to produce the 
observed dust properties in these quasars. This suggests that 
the quasar hosts should have already built up a large stellar 
population by z ~ 6. Our results support a merger-driven ori- 
gin for luminous quasars at high redshifts, and provide further 
evidence for the hypothesis of starburst-to-quasar evolution. 
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